v0.14.0
Loading...
Searching...
No Matches
EdgeForce.cpp
Go to the documentation of this file.
1
2/** \file EdgeForce.cpp
3 \ingroup mofem_static_boundary_conditions
4*/
5
6
7
8#include <MoFEM.hpp>
9using namespace MoFEM;
11#include <SurfacePressure.hpp>
12#include <EdgeForce.hpp>
13
15 const std::string field_name, Vec f, bCForce &data,
16 boost::ptr_vector<MethodForForceScaling> &methods_op, bool use_snes_f)
18 F(f), dAta(data), methodsOp(methods_op), useSnesF(use_snes_f) {}
19
24
25 if (data.getIndices().size() == 0) {
27 }
28 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
29 if (dAta.eDges.find(ent) == dAta.eDges.end()) {
31 }
32
33 // Get pointer to DOF and its rank
34 const auto &dof_ptr = data.getFieldDofs()[0];
35 int rank = dof_ptr->getNbOfCoeffs();
36
37 int nb_dofs = data.getIndices().size();
38
39 Nf.resize(nb_dofs, false);
40 Nf.clear();
41
42 int nb_gauss_pts = data.getN().size1();
43 wEights.resize(nb_gauss_pts, false);
44
45 // This will work for fluxes and other fields with rank other than 3.
46 for (int rr = 0; rr < rank; rr++) {
47
48 // Get force value for each vector element from blockset data.
49 double force;
50 if (rr == 0) {
51 force = dAta.data.data.value3 * dAta.data.data.value1;
52 } else if (rr == 1) {
53 force = dAta.data.data.value4 * dAta.data.data.value1;
54 } else if (rr == 2) {
55 force = dAta.data.data.value5 * dAta.data.data.value1;
56 } else {
57 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
58 }
59
60 // Integrate force on the line
61 for (int gg = 0; gg < nb_gauss_pts; gg++) {
62
63 if (!rr) {
64 wEights[gg] = 0;
65 // This is if edge is curved, i.e. HO geometry
66 for (int dd = 0; dd < 3; dd++) {
67 wEights[gg] += pow(getTangentAtGaussPts()(gg, dd), 2);
68 }
69 wEights[gg] = std::sqrt(wEights[gg]);
70 wEights[gg] *= getGaussPts()(1, gg);
71 }
72
73 cblas_daxpy(nb_dofs / rank, wEights[gg] * force, &data.getN()(gg, 0), 1,
74 &Nf[rr], rank);
75 }
76 }
77
78 // I time/step varying force or calculate in arc-length control. This hack
79 // scale force appropriately, and is controlled for user
81
82 // Assemble force into right-hand vector
83 Vec myF = F;
84 if (useSnesF || F == PETSC_NULL)
85 myF = getKSPf();
86
87 CHKERR VecSetValues(myF, data.getIndices().size(), &data.getIndices()[0],
88 &Nf[0], ADD_VALUES);
89
91}
92
94 int ms_id, bool use_snes_f) {
95 const CubitMeshSets *cubit_meshset_ptr;
96 MeshsetsManager *mmanager_ptr;
98 CHKERR mField.getInterface(mmanager_ptr);
99 CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, NODESET, &cubit_meshset_ptr);
100 CHKERR cubit_meshset_ptr->getBcDataStructure(mapForce[ms_id].data);
101 CHKERR mField.get_moab().get_entities_by_type(
102 cubit_meshset_ptr->meshset, MBEDGE, mapForce[ms_id].eDges, true);
103 // Add operator for element, set data and entities operating on the data
104 fe.getOpPtrVector().push_back(
105 new OpEdgeForce(field_name, F, mapForce[ms_id], methodsOp, use_snes_f));
107}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#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
@ NODESET
Definition: definitions.h:146
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: EdgeForce.cpp:21
OpEdgeForce(const std::string field_name, Vec f, bCForce &data, boost::ptr_vector< MethodForForceScaling > &methods_op, bool use_snes_f=false)
Definition: EdgeForce.cpp:14
MoFEM::Interface & mField
Definition: EdgeForce.hpp:15
MoFEMErrorCode addForce(const std::string field_name, Vec F, int ms_id, bool use_snes_f=false)
Definition: EdgeForce.cpp:93
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: EdgeForce.hpp:34
std::map< int, bCForce > mapForce
Definition: EdgeForce.hpp:32
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
this struct keeps basic methods for moab meshset about material and boundary conditions
MoFEMErrorCode getBcDataStructure(CUBIT_BC_DATA_TYPE &data) const
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.