v0.9.1
BodyForce.hpp
Go to the documentation of this file.
1 /** \file BodyForce.hpp
2  */
3 
4 /* This file is part of MoFEM.
5  * MoFEM is free software: you can redistribute it and/or modify it under
6  * the terms of the GNU Lesser General Public License as published by the
7  * Free Software Foundation, either version 3 of the License, or (at your
8  * option) any later version.
9  *
10  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13  * License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
17 
18 #ifndef __BODY_FORCE_HPP
19 #define __BODY_FORCE_HPP
20 
21 /** \brief Body forces elements
22  * \ingroup mofem_body_forces
23  */
25 
27 
31  int getRule(int order) { return order; };
32  };
33 
35  MyVolumeFE &getLoopFe() { return fe; }
36 
38  : mField(m_field), fe(m_field) {}
39 
40  struct OpBodyForce
42 
43  Vec F;
44  Block_BodyForces &dAta;
45  Range blockTets;
46  OpBodyForce(const std::string field_name, Vec _F, Block_BodyForces &data,
47  Range block_tets)
49  field_name, UserDataOperator::OPROW),
50  F(_F), dAta(data), blockTets(block_tets) {}
51 
53 
54  MoFEMErrorCode doWork(int side, EntityType type,
57 
58  if (data.getIndices().size() == 0)
60  if (blockTets.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
61  blockTets.end())
63 
64  const auto &dof_ptr = data.getFieldDofs()[0];
65  int rank = dof_ptr->getNbOfCoeffs();
66  int nb_row_dofs = data.getIndices().size() / rank;
67 
68  Nf.resize(data.getIndices().size());
69  bzero(&*Nf.data().begin(), data.getIndices().size() * sizeof(FieldData));
70 
71  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
72  double val = getVolume() * getGaussPts()(3, gg);
73  if (getHoGaussPtsDetJac().size() > 0) {
74  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
75  }
76  for (int rr = 0; rr < rank; rr++) {
77 
78  double acc;
79  if (rr == 0) {
80  acc = -dAta.data.acceleration_x;
81  } else if (rr == 1) {
82  acc = -dAta.data.acceleration_y;
83  } else if (rr == 2) {
84  acc = -dAta.data.acceleration_z;
85  } else {
86  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
87  }
88  acc *= dAta.data.density;
89  cblas_daxpy(nb_row_dofs, val * acc, &data.getN()(gg, 0), 1, &Nf[rr],
90  rank);
91  }
92  }
93 
94  CHKERR VecSetValues(F, data.getIndices().size(), &data.getIndices()[0],
95  &Nf[0], ADD_VALUES);
96 
98  }
99  };
100 
101  MoFEMErrorCode addBlock(const std::string field_name, Vec &F, int ms_id) {
102  const CubitMeshSets *cubit_meshset_ptr;
103  MeshsetsManager *mmanager_ptr;
105  CHKERR mField.getInterface(mmanager_ptr);
106  CHKERR mmanager_ptr->getCubitMeshsetPtr(ms_id, BLOCKSET,
107  &cubit_meshset_ptr);
108  CHKERR cubit_meshset_ptr->getAttributeDataStructure(mapData[ms_id]);
109  EntityHandle meshset = cubit_meshset_ptr->getMeshset();
110  Range tets;
111  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets, true);
112  fe.getOpPtrVector().push_back(
113  new OpBodyForce(field_name, F, mapData[ms_id], tets));
115  }
116 
117 private:
118  std::map<int, Block_BodyForces> mapData;
119 };
120 
121 /// \brief USe BodyForceConstantField
123 
124 #endif //__BODY_FORCE_HPP
125 
126 /**
127  * \defgroup mofem_body_forces Body forces elements
128  * \ingroup user_modules
129  */
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Deprecated interface functions.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
virtual moab::Interface & get_moab()=0
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
std::map< int, Block_BodyForces > mapData
Definition: BodyForce.hpp:118
MoFEMErrorCode addBlock(const std::string field_name, Vec &F, int ms_id)
Definition: BodyForce.hpp:101
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
DEPRECATED typedef BodyForceConstantField BodyFroceConstantField
USe BodyForceConstantField.
Definition: BodyForce.hpp:122
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
MyVolumeFE(MoFEM::Interface &m_field)
Definition: BodyForce.hpp:29
MoFEM::Interface & mField
Definition: BodyForce.hpp:26
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Body forces elements.
Definition: BodyForce.hpp:24
constexpr int order
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: BodyForce.hpp:54
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:601
double FieldData
Field data type.
Definition: Types.hpp:36
MyVolumeFE & getLoopFe()
Definition: BodyForce.hpp:35
#define DEPRECATED
Definition: definitions.h:29
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
#define _F(n)
Definition: quad.c:25
OpBodyForce(const std::string field_name, Vec _F, Block_BodyForces &data, Range block_tets)
Definition: BodyForce.hpp:46
BodyForceConstantField(MoFEM::Interface &m_field)
Definition: BodyForce.hpp:37