33 #include "libmesh/parsed_function.h"
34 #include "libmesh/zero_function.h"
44 _quadratic_scaling(false)
60 std::string base_physics_name =
"VelocityPenalty";
63 base_physics_name.push_back(
'2');
67 base_physics_name.push_back(
'3');
69 std::string penalty_function =
70 input(
"Physics/"+base_physics_name+
"/penalty_function",
73 if (penalty_function ==
"0")
74 this->normal_vector_function.reset
75 (
new libMesh::ZeroFunction<libMesh::Number>());
77 this->normal_vector_function.reset
78 (
new libMesh::ParsedFunction<libMesh::Number>(penalty_function));
80 std::string base_function =
81 input(
"Physics/"+base_physics_name+
"/base_velocity",
84 if (penalty_function ==
"0" && base_function ==
"0")
85 std::cout <<
"Warning! Zero VelocityPenalty specified!" << std::endl;
87 if (base_function ==
"0")
88 this->base_velocity_function.reset
89 (
new libMesh::ZeroFunction<libMesh::Number>());
91 this->base_velocity_function.reset
92 (
new libMesh::ParsedFunction<libMesh::Number>(base_function));
95 input(
"Physics/"+base_physics_name+
"/quadratic_scaling",
false);
100 (
const libMesh::Point& point,
101 const libMesh::Real time,
102 const libMesh::NumberVectorValue& U,
103 libMesh::NumberVectorValue& F,
104 libMesh::NumberTensorValue *dFdU)
109 libmesh_assert(normal_vector_function.get());
110 libmesh_assert(base_velocity_function.get());
112 libMesh::DenseVector<libMesh::Number> output_vec(3);
114 (*normal_vector_function)(point, time,
117 libMesh::NumberVectorValue U_N(output_vec(0),
121 (*base_velocity_function)(point, time,
124 const libMesh::NumberVectorValue U_B(output_vec(0),
128 const libMesh::NumberVectorValue U_Rel = U-U_B;
134 const libMesh::Number U_N_mag = std::sqrt(U_N*U_N);
139 const libMesh::NumberVectorValue U_N_unit = U_N/U_N_mag;
141 F = -(U_Rel*U_N)*U_N_unit;
144 for (
unsigned int i=0; i != 3; ++i)
145 for (
unsigned int j=0; j != 3; ++j)
146 (*dFdU)(i,j) = -(U_N(j))*U_N_unit(i);
149 if (_quadratic_scaling)
151 const libMesh::Number U_Rel_mag = std::sqrt(U_Rel * U_Rel);
160 const libMesh::NumberVectorValue U_Rel_unit = U_Rel/U_Rel_mag;
162 (*dFdU) *= U_Rel_mag;
165 for (
unsigned int i=0; i != 3; ++i)
166 for (
unsigned int j=0; j != 3; ++j)
167 (*dFdU)(i,j) += F(i)*U_Rel_unit(j);
174 bool do_correction =
false;
178 libmesh_not_implemented();
180 const libMesh::Number U_Rel_mag_sq = U_Rel * U_Rel;
183 F -= (U_Rel*F)*U_Rel/U_Rel_mag_sq;
const PhysicsName incompressible_navier_stokes
Physics class for Incompressible Navier-Stokes.
const PhysicsName velocity_penalty2_adjoint_stab
INSTANTIATE_INC_NS_SUBCLASS(VelocityPenaltyBase)
const PhysicsName velocity_penalty3
const PhysicsName velocity_penalty2
virtual void read_input_options(const GetPot &input)
Read options from GetPot input file.
const PhysicsName velocity_penalty3_adjoint_stab
bool compute_force(const libMesh::Point &point, const libMesh::Real time, const libMesh::NumberVectorValue &U, libMesh::NumberVectorValue &F, libMesh::NumberTensorValue *dFdU=NULL)