v0.14.0
Loading...
Searching...
No Matches
FluidPressure.cpp
Go to the documentation of this file.
1/* \file FluidPressure.cpp
2 *
3 * \brief Implementation of fluid pressure element
4 *
5 * \todo Implement nonlinear case (consrvative force, i.e. normal follows
6 * surface normal)
7 *
8 */
9
10
11
13 int side, EntityType type, EntitiesFieldData::EntData &data) {
15 if (data.getIndices().size() == 0)
18 if (dAta.tRis.find(ent) == dAta.tRis.end())
20
21 const auto &dof_ptr = data.getFieldDofs()[0];
22 const int rank = dof_ptr->getNbOfCoeffs();
23 const int nb_row_dofs = data.getIndices().size() / rank;
24
25 Nf.resize(data.getIndices().size());
26 Nf.clear();
27
28 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
29
30 VectorDouble dist;
31 VectorDouble zero_pressure = dAta.zEroPressure;
32
33 dist = ublas::matrix_row<MatrixDouble>(getCoordsAtGaussPts(), gg);
34 dist -= zero_pressure;
35 double dot = cblas_ddot(3, &dist[0], 1, &dAta.aCCeleration[0], 1);
37 dot = fmax(0, dot);
38 double pressure = dot * dAta.dEnsity;
39
40 for (int rr = 0; rr < rank; rr++) {
41 double force;
42 if (hoGeometry) {
43 force = pressure * getNormalsAtGaussPts()(gg, rr);
44 } else {
45 force = pressure * getNormal()[rr];
46 }
47 cblas_daxpy(nb_row_dofs, getGaussPts()(2, gg) * force,
48 &data.getN()(gg, 0), 1, &Nf[rr], rank);
49 }
50 }
51
52 // Scale force using user defined scaling operator
54
55 auto get_f = [&]() {
56 if (F == PETSC_NULL)
57 return getKSPf();
58 return F;
59 };
60
61 // Assemble force into vector
62 CHKERR VecSetValues(get_f(), data, &*Nf.data().begin(), ADD_VALUES);
63
65}
66
68 const std::string field_name, const std::string mesh_nodals_positions) {
70
71 CHKERR mField.add_finite_element("FLUID_PRESSURE_FE", MF_ZERO);
78 if (mField.check_field(mesh_nodals_positions)) {
80 mesh_nodals_positions);
81 }
82
83 // takes skin of block of entities
84 Skinner skin(&mField.get_moab());
85 // loop over all blocksets and get data which name is FluidPressure
87
88 if (bit->getName().compare(0, 14, "FLUID_PRESSURE") == 0) {
89
90 // get block attributes
91 std::vector<double> attributes;
92 CHKERR bit->getAttributes(attributes);
93 if (attributes.size() < 7) {
94 SETERRQ1(PETSC_COMM_SELF, 1,
95 "not enough block attributes to deffine fluid pressure "
96 "element, attributes.size() = %d ",
97 attributes.size());
98 }
99 setOfFluids[bit->getMeshsetId()].dEnsity = attributes[0];
100 setOfFluids[bit->getMeshsetId()].aCCeleration.resize(3);
101 setOfFluids[bit->getMeshsetId()].aCCeleration[0] = attributes[1];
102 setOfFluids[bit->getMeshsetId()].aCCeleration[1] = attributes[2];
103 setOfFluids[bit->getMeshsetId()].aCCeleration[2] = attributes[3];
104 setOfFluids[bit->getMeshsetId()].zEroPressure.resize(3);
105 setOfFluids[bit->getMeshsetId()].zEroPressure[0] = attributes[4];
106 setOfFluids[bit->getMeshsetId()].zEroPressure[1] = attributes[5];
107 setOfFluids[bit->getMeshsetId()].zEroPressure[2] = attributes[6];
108 // get blok tetrahedron and triangles
109 Range tets;
110 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTET, tets,
111 true);
112 Range tris;
113 CHKERR mField.get_moab().get_entities_by_type(
114 bit->meshset, MBTRI, setOfFluids[bit->getMeshsetId()].tRis, true);
115 // this get triangles only on block surfaces
116 Range tets_skin_tris;
117 CHKERR skin.find_skin(0, tets, false, tets_skin_tris);
118 setOfFluids[bit->getMeshsetId()].tRis.merge(tets_skin_tris);
119 std::ostringstream ss;
120 ss << setOfFluids[bit->getMeshsetId()] << std::endl;
121 PetscPrintf(mField.get_comm(), ss.str().c_str());
122
124 setOfFluids[bit->getMeshsetId()].tRis, MBTRI, "FLUID_PRESSURE_FE");
125 }
126 }
127
129}
130
132 string field_name, Vec F, bool allow_negative_pressure, bool ho_geometry) {
134 std::map<MeshSetId, FluidData>::iterator sit = setOfFluids.begin();
135 for (; sit != setOfFluids.end(); sit++) {
136 // add finite element
137 fe.getOpPtrVector().push_back(
138 new OpCalculatePressure(field_name, F, sit->second, methodsOp,
139 allow_negative_pressure, ho_geometry));
140 }
142}
143
144std::ostream &operator<<(std::ostream &os, const FluidPressure::FluidData &e) {
145 os << "dEnsity " << e.dEnsity << std::endl;
146 os << "aCCeleration " << e.aCCeleration << std::endl;
147 os << "zEroPressure " << e.zEroPressure << std::endl;
148 return os;
149}
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
@ F
auto get_f
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
constexpr auto field_name
VectorDouble zEroPressure
fluid level of reference zero pressure.
Range tRis
range of surface element to which fluid pressure is applied
double dEnsity
fluid density [kg/m^2] or any consistent unit
VectorDouble aCCeleration
acceleration [m/s^2]
boost::ptr_vector< MethodForForceScaling > & methodsOp
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
bool allowNegativePressure
allows for negative pressures
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< MeshSetId, FluidData > setOfFluids
MoFEMErrorCode setNeumannFluidPressureFiniteElementOperators(string field_name, Vec F, bool allow_negative_pressure=true, bool ho_geometry=false)
MyTriangleFE fe
MoFEM::Interface & mField
MoFEMErrorCode addNeumannFluidPressureBCElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.