v0.13.1
Loading...
Searching...
No Matches
DirichletBC.hpp
Go to the documentation of this file.
1/* \file Dirichlet.hpp
2 * \brief Implementation of Dirichlet boundary conditions
3 * \ingroup Dirichlet_bc
4 *
5 * Structures and method in this file erase rows and column, set value on
6 * matrix diagonal and on the right hand side vector to enforce boundary
7 * condition.
8 *
9 * Current implementation is suboptimal, classes name too long. Need to
10 * rethinking and improved, more elegant and more efficient implementation.
11 *
12 */
13
14/* Notes:
15
16 DirichletSetFieldFromBlock implemented by Zahur Ullah
17 (Zahur.Ullah@glasgow.ac.uk)
18
19 */
20
21
22
23#ifndef __DIRICHLET_HPP__
24#define __DIRICHLET_HPP__
25
26using namespace boost::numeric;
27
28/** \brief Data from Cubit blocksets
29 * \ingroup Dirichlet_bc
30 */
31struct DataFromBc {
32 VectorDouble scaled_values;
33 VectorDouble initial_values;
34 VectorInt bc_flags;
36
37 // for rotation
41 double theta;
42
44 : scaled_values(3), initial_values(3), bc_flags(3), is_rotation(false) {}
45
46 MoFEMErrorCode getBcData(DisplacementCubitBcData &mydata,
47 const MoFEM::CubitMeshSets *it);
48 MoFEMErrorCode getBcData(std::vector<double> &mydata,
49 const MoFEM::CubitMeshSets *it);
50 MoFEMErrorCode getEntitiesFromBc(MoFEM::Interface &mField,
51 const MoFEM::CubitMeshSets *it);
52};
53
54/** \brief Set Dirichlet boundary conditions on displacements
55 * \ingroup Dirichlet_bc
56 */
58
60 const std::string fieldName; ///< field name to set Dirichlet BC
61 double dIag; ///< diagonal value set on zeroed column and rows
62
64 const std::string &field_name, Mat Aij, Vec X, Vec F,
65 string blockset_name = "DISPLACEMENT");
67 const std::string &field_name,
68 string blockset_name = "DISPLACEMENT");
69
70 std::map<DofIdx, FieldData> mapZeroRows;
71 std::vector<int> dofsIndices;
72 std::vector<double> dofsValues;
73 std::vector<double> dofsXValues;
74 const std::string blocksetName;
75
76 boost::ptr_vector<MethodForForceScaling> methodsOp;
77 virtual MoFEMErrorCode iNitialize();
78
79 MoFEMErrorCode preProcess();
80 MoFEMErrorCode operator()() { return 0; }
81 MoFEMErrorCode postProcess();
82 /**
83 * @brief Get the Bc Data From Sets And Blocks object
84 * Use DISPLACEMENT blockset name (default)
85 * with 6 atributes:
86 * 1,2,3 are values of displacements x,y,z
87 * 4,5,6 are flags for x,y,z (0 or 1)
88 * @param bc_data
89 * @return MoFEMErrorCode
90 */
91 MoFEMErrorCode getBcDataFromSetsAndBlocks(std::vector<DataFromBc> &bc_data);
92 /**
93 * @brief Get the Rotation Bc From Block object
94 * Use ROTATION blockset name
95 * with 7 atributes:
96 * 1,2,3 are x,y,z coords of the center of rotation
97 * 4,5,6 are are angular velocities in x,y,z
98 * @param bc_data
99 * @return MoFEMErrorCode
100 */
101 MoFEMErrorCode getRotationBcFromBlock(std::vector<DataFromBc> &bc_data);
102
103 /**
104 * @brief Calculate displacements from rotation for particular dof
105 * @param dof
106 * @param bc_data
107 * @return MoFEMErrorCode
108 */
109 MoFEMErrorCode calculateRotationForDof(VectorDouble3 &coords,
110 DataFromBc &bc_data);
111 MoFEMErrorCode calculateRotationForDof(EntityHandle ent,
112 DataFromBc &bc_data);
113 MoFEMErrorCode applyScaleBcData(DataFromBc &bc_data);
114};
115
120 DataFromBc &data_from_dirichlet_bc)
121 : dirichletBcPtr(dirichlet_bc_ptr), dataFromDirichletBc(data_from_dirichlet_bc) {}
122
123 MoFEMErrorCode preProcess() { return 0; }
124 MoFEMErrorCode postProcess() { return 0; }
125 MoFEMErrorCode operator()() {
127 auto &mField = dirichletBcPtr->mField;
128 auto &bc_it = dataFromDirichletBc;
129
130 EntityHandle v = entPtr->getEnt();
131 int coeff = fieldPtr->getNbOfCoeffs();
133 for (int i = 0; i != coeff; i++) {
134 if (bc_it.bc_flags[i]) {
135 if (entPtr->getEntType() == MBVERTEX) {
136 entPtr->getEntFieldData()[i] = bc_it.scaled_values[i];
137 } else if (!entPtr->getEntFieldData().empty()) {
138 entPtr->getEntFieldData()[i] = 0;
139 }
140 }
141 }
142
144 }
145};
146
148 // using BcEntMethodDisp::BcEntMethodDisp;
153 DataFromBc &data_from_dirichlet_bc,
154 string material_positions)
155 : dirichletBcPtr(dirichlet_bc_ptr),
156 dataFromDirichletBc(data_from_dirichlet_bc),
157 materialPositions(material_positions) {}
158
159 MoFEMErrorCode preProcess() { return 0; }
160 MoFEMErrorCode postProcess() { return 0; }
161
162 MoFEMErrorCode operator()() {
164 EntityHandle ent = entPtr->getEnt();
165 auto &mField = dirichletBcPtr->mField;
166 auto &bc_it = dataFromDirichletBc;
167 EntityHandle v = entPtr->getEnt();
168
169 const FieldEntity_multiIndex *field_ents;
170 CHKERR mField.get_field_ents(&field_ents);
171 auto &field_ents_by_uid = field_ents->get<Unique_mi_tag>();
172
173 auto get_coords = [&]() {
174 VectorDouble coords({0, 0, 0});
175 if (entPtr->getEntType() == MBVERTEX) {
176 auto eit =
177 field_ents_by_uid.find(FieldEntity::getLocalUniqueIdCalculate(
178 mField.get_field_bit_number(materialPositions), ent));
179 if (eit != field_ents_by_uid.end())
180 noalias(coords) = (*eit)->getEntFieldData();
181 else
182 CHKERR mField.get_moab().get_coords(&ent, 1, &*coords.data().begin());
183 }
184 return coords;
185 };
186
187 int coeff = fieldPtr->getNbOfCoeffs();
188 auto coords = get_coords();
189
191 for (int i = 0; i != coeff; i++) {
192 if (bc_it.bc_flags[i]) {
193 if (entPtr->getEntType() == MBVERTEX) {
194 entPtr->getEntFieldData()[i] = coords(i) + bc_it.scaled_values[i];
195 } else if (!entPtr->getEntFieldData().empty()) {
196 entPtr->getEntFieldData()[i] = 0;
197 }
198 }
199 }
200
202 }
203};
204
205/// \deprecated use DirichletDisplacementBc
207
208/** \brief Set Dirichlet boundary conditions on spatial displacements
209 * \ingroup Dirichlet_bc
210 */
212
214 MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x,
215 Vec f, const std::string material_positions = "MESH_NODE_POSITIONS",
216 const std::string blockset_name = "DISPLACEMENT")
217 : DirichletDisplacementBc(m_field, field_name, aij, x, f, blockset_name),
218 materialPositions(material_positions) {}
219
221 MoFEM::Interface &m_field, const std::string &field_name,
222 const std::string material_positions = "MESH_NODE_POSITIONS",
223 const std::string blockset_name = "DISPLACEMENT")
224 : DirichletDisplacementBc(m_field, field_name, blockset_name),
225 materialPositions(material_positions) {}
226
227 std::string materialPositions; ///< name of the field with reference material
228 ///< positions
229 std::vector<std::string> fixFields; ///<
230
231 VectorDouble cOords;
232 MoFEMErrorCode iNitialize();
233};
234
235/// \deprecated use DirichletSpatialPositionsBc
238
240
242 const std::string &field_name, Mat aij, Vec x, Vec f)
243 : DirichletDisplacementBc(m_field, field_name, aij, x, f) {}
244
246 const std::string &field_name)
248
249 MoFEMErrorCode iNitialize();
250};
251
252/// \deprecated use DirichletTemperatureBc
254
255/** \brief Fix dofs on entities
256 * \ingroup Dirichlet_bc
257 */
259
261 std::vector<std::string> fieldNames;
263 const std::string field_name, Mat aij, Vec x,
264 Vec f, Range &ents)
265 : DirichletDisplacementBc(m_field, field_name, aij, x, f), eNts(ents) {
266 fieldNames.push_back(fieldName);
267 }
268
270 const std::string field_name, Range &ents)
271 : DirichletDisplacementBc(m_field, field_name), eNts(ents) {
272 fieldNames.push_back(fieldName);
273 }
274
275 MoFEMErrorCode iNitialize();
276 MoFEMErrorCode preProcess();
277 MoFEMErrorCode postProcess();
278};
279
280/** \brief Set Dirichlet boundary conditions on displacements by removing dofs
281 * \ingroup Dirichlet_bc
282 */
284
285 boost::shared_ptr<vector<DataFromBc>> bcDataPtr;
288
290 const std::string &field_name,
291 const std::string &problem_name,
292 string blockset_name = "DISPLACEMENT",
293 bool is_partitioned = false)
294 : DirichletDisplacementBc(m_field, field_name, blockset_name),
295 problemName(problem_name), isPartitioned(is_partitioned) {}
296
297 MoFEMErrorCode iNitialize();
298
299 virtual boost::shared_ptr<EntityMethod> getEntMethodPtr(DataFromBc &data) {
300 return boost::make_shared<BcEntMethodDisp>(this, data);
301 }
302
303 MoFEMErrorCode preProcess();
304 MoFEMErrorCode operator()() { return 0; }
305 MoFEMErrorCode postProcess() { return 0; }
306};
307
308/** \brief Set Dirichlet boundary conditions on spatial positions by removing dofs
309 * \ingroup Dirichlet_bc
310 */
312
313 std::string materialPositions;
314
316 MoFEM::Interface &m_field, const std::string &field_name,
317 const std::string &problem_name,
318 const std::string material_positions = "MESH_NODE_POSITIONS",
319 string blockset_name = "DISPLACEMENT", bool is_partitioned = false)
320 : DirichletDisplacementRemoveDofsBc(m_field, field_name, problem_name, blockset_name,
321 is_partitioned),
322 materialPositions(material_positions) {}
323
324 boost::shared_ptr<EntityMethod> getEntMethodPtr(DataFromBc &data) override {
325 return boost::make_shared<BcEntMethodSpatial>(this, data,
327 }
328};
329
330/// \deprecated use DirichletFixFieldAtEntitiesBc
332
333/**
334 * \brief Add boundary conditions form block set having 6 attributes
335 *
336 * First 3 values are magnitudes of dofs e.g. in x,y,z direction and next 3 are
337 flags, respectively.
338 * If flag is false ( = 0), particular dof is not taken into account.
339 Usage in Cubit for displacement:
340 block 1 tri 28 32
341 block 1 name "DISPLACEMENT_1"
342 block 1 attribute count 6
343 block 1 attribute index 1 97 # any value
344 block 1 attribute index 2 0
345 block 1 attribute index 3 0
346 block 1 attribute index 4 0 # flag for x dir
347 block 1 attribute index 5 1 # flag for y dir
348 block 1 attribute index 6 1 # flag for z dir
349 This means that we set zero displacement on y and z dir and on x set
350 direction freely. (value 97 is irrelevant because flag for 1 value is 0
351 (false)) It can be usefull if we want to set boundary conditions directly to
352 triangles e.g, since standard boundary conditions in Cubit allow only using
353 nodeset or surface which might not work with mesh based on facet engine (e.g.
354 STL file)
355 */
357
359 const std::string &field_name,
360 const std::string &blockset_name, Mat aij,
361 Vec x, Vec f)
362 : DirichletDisplacementBc(m_field, field_name, aij, x, f, blockset_name) {
363 }
364
366 const std::string &field_name,
367 const std::string &blockset_name)
368 : DirichletDisplacementBc(m_field, field_name, blockset_name) {}
369};
370
371/// \deprecated use DirichletSetFieldFromBlockWithFlags
374
375/// \deprecated use DirichletSetFieldFromBlockWithFlags
378
379/// \deprecated use DirichletSetFieldFromBlockWithFlags
382/**
383 * @brief calculate reactions from vector of internal forces on meshsets
384 *
385 * example usage
386 *
387 * \code
388 Vec F_int;
389 DMCreateGlobalVector_MoFEM(dm, &F_int);
390
391 feRhs->snes_ctx = FEMethod::CTX_SNESSETFUNCTION;
392 feRhs->snes_f = F_int;
393 DMoFEMLoopFiniteElements(dm, "ELASTIC", feRhs);
394
395 VecAssemblyBegin(F_int);
396 VecAssemblyEnd(F_int);
397 VecGhostUpdateBegin(F_int, INSERT_VALUES, SCATTER_FORWARD);
398 VecGhostUpdateEnd(F_int, INSERT_VALUES, SCATTER_FORWARD);
399
400 Reactions my_react(m_field, "DM_ELASTIC", "U");
401 my_react.calculateReactions(F_int);
402 int fix_nodes_meshset_id = 1;
403 cout << my_react.getReactionsFromSet(fix_nodes_meshset_id) << endl;
404
405* \endcode
406 */
407struct Reactions {
408
409 Reactions(MoFEM::Interface &m_field, string problem_name, string field_name)
410 : mField(m_field), problemName(problem_name), fieldName(field_name) {}
411
412 typedef std::map<int, VectorDouble> ReactionsMap;
414 /**
415 * @brief Get the Reactions Map
416 *
417 * @return const ReactionsMap&
418 */
419 inline const ReactionsMap &getReactionsMap() const { return reactionsMap; }
420 /**
421 * @brief Get the Reactions at specified meshset id
422 *
423 * @param id meshset id (from Cubit)
424 * @return const VectorDouble&
425 */
426 inline const VectorDouble &getReactionsFromSet(const int &id) const {
427 return reactionsMap.at(id);
428 }
429 /**
430 * @brief calculate reactions from a given vector
431 *
432 * @param internal forces vector
433 * @return MoFEMErrorCode
434 */
435 MoFEMErrorCode calculateReactions(Vec &internal);
436
437private:
438 std::string problemName;
439 std::string fieldName;
441};
442
443#endif //__DIRICHLET_HPP__
444
445/**
446 * \defgroup Dirichlet_bc Dirichlet boundary conditions
447 * \ingroup user_modules
448 **/
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletSetFieldFromBlock
DEPRECATED typedef DirichletFixFieldAtEntitiesBc FixBcAtEntities
DEPRECATED typedef DirichletDisplacementBc DisplacementBCFEMethodPreAndPostProc
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletBCFromBlockSetFEMethodPreAndPostProc
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletBCFromBlockSetFEMethodPreAndPostProcWithFlags
DEPRECATED typedef DirichletTemperatureBc TemperatureBCFEMethodPreAndPostProc
DEPRECATED typedef DirichletSpatialPositionsBc SpatialPositionsBCFEMethodPreAndPostProc
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define DEPRECATED
Definition: definitions.h:17
#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
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, member< FieldEntity, UId, &FieldEntity::localUId > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity::interface_type_RefEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex
MultiIndex container keeps FieldEntity.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
constexpr auto field_name
MoFEMErrorCode postProcess()
function is run at the end of loop
DirichletDisplacementBc * dirichletBcPtr
MoFEMErrorCode operator()()
function is run for every finite element
BcEntMethodDisp(DirichletDisplacementBc *dirichlet_bc_ptr, DataFromBc &data_from_dirichlet_bc)
DataFromBc & dataFromDirichletBc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
DataFromBc & dataFromDirichletBc
MoFEMErrorCode postProcess()
function is run at the end of loop
BcEntMethodSpatial(DirichletDisplacementBc *dirichlet_bc_ptr, DataFromBc &data_from_dirichlet_bc, string material_positions)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
DirichletDisplacementBc * dirichletBcPtr
MoFEMErrorCode operator()()
function is run for every finite element
Data from Cubit blocksets.
Definition: DirichletBC.hpp:31
VectorInt bc_flags
Definition: DirichletBC.hpp:34
FTensor::Tensor1< double, 3 > t_centr
Definition: DirichletBC.hpp:40
MoFEMErrorCode getBcData(DisplacementCubitBcData &mydata, const MoFEM::CubitMeshSets *it)
MoFEMErrorCode getEntitiesFromBc(MoFEM::Interface &mField, const MoFEM::CubitMeshSets *it)
VectorDouble scaled_values
Definition: DirichletBC.hpp:32
FTensor::Tensor1< double, 3 > t_normal
Definition: DirichletBC.hpp:39
VectorDouble initial_values
Definition: DirichletBC.hpp:33
double theta
Definition: DirichletBC.hpp:41
bool is_rotation
Definition: DirichletBC.hpp:38
Range bc_ents[3]
Definition: DirichletBC.hpp:35
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
const std::string blocksetName
Definition: DirichletBC.hpp:74
MoFEMErrorCode postProcess()
function is run at the end of loop
std::map< DofIdx, FieldData > mapZeroRows
Definition: DirichletBC.hpp:70
MoFEMErrorCode getRotationBcFromBlock(std::vector< DataFromBc > &bc_data)
Get the Rotation Bc From Block object Use ROTATION blockset name with 7 atributes: 1,...
const std::string fieldName
field name to set Dirichlet BC
Definition: DirichletBC.hpp:60
MoFEMErrorCode operator()()
function is run for every finite element
Definition: DirichletBC.hpp:80
double dIag
diagonal value set on zeroed column and rows
Definition: DirichletBC.hpp:61
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:71
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode getBcDataFromSetsAndBlocks(std::vector< DataFromBc > &bc_data)
Get the Bc Data From Sets And Blocks object Use DISPLACEMENT blockset name (default) with 6 atributes...
Definition: DirichletBC.cpp:64
MoFEMErrorCode applyScaleBcData(DataFromBc &bc_data)
std::vector< double > dofsXValues
Definition: DirichletBC.hpp:73
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: DirichletBC.hpp:76
std::vector< double > dofsValues
Definition: DirichletBC.hpp:72
MoFEM::Interface & mField
Definition: DirichletBC.hpp:59
virtual MoFEMErrorCode iNitialize()
MoFEMErrorCode calculateRotationForDof(VectorDouble3 &coords, DataFromBc &bc_data)
Calculate displacements from rotation for particular dof.
Set Dirichlet boundary conditions on displacements by removing dofs.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual boost::shared_ptr< EntityMethod > getEntMethodPtr(DataFromBc &data)
DirichletDisplacementRemoveDofsBc(MoFEM::Interface &m_field, const std::string &field_name, const std::string &problem_name, string blockset_name="DISPLACEMENT", bool is_partitioned=false)
boost::shared_ptr< vector< DataFromBc > > bcDataPtr
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode postProcess()
function is run at the end of loop
Fix dofs on entities.
std::vector< std::string > fieldNames
DirichletFixFieldAtEntitiesBc(MoFEM::Interface &m_field, const std::string field_name, Mat aij, Vec x, Vec f, Range &ents)
DirichletFixFieldAtEntitiesBc(MoFEM::Interface &m_field, const std::string field_name, Range &ents)
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Add boundary conditions form block set having 6 attributes.
DirichletSetFieldFromBlockWithFlags(MoFEM::Interface &m_field, const std::string &field_name, const std::string &blockset_name)
DirichletSetFieldFromBlockWithFlags(MoFEM::Interface &m_field, const std::string &field_name, const std::string &blockset_name, Mat aij, Vec x, Vec f)
Set Dirichlet boundary conditions on spatial displacements.
MoFEMErrorCode iNitialize()
std::vector< std::string > fixFields
DirichletSpatialPositionsBc(MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x, Vec f, const std::string material_positions="MESH_NODE_POSITIONS", const std::string blockset_name="DISPLACEMENT")
DirichletSpatialPositionsBc(MoFEM::Interface &m_field, const std::string &field_name, const std::string material_positions="MESH_NODE_POSITIONS", const std::string blockset_name="DISPLACEMENT")
Set Dirichlet boundary conditions on spatial positions by removing dofs.
boost::shared_ptr< EntityMethod > getEntMethodPtr(DataFromBc &data) override
DirichletSpatialRemoveDofsBc(MoFEM::Interface &m_field, const std::string &field_name, const std::string &problem_name, const std::string material_positions="MESH_NODE_POSITIONS", string blockset_name="DISPLACEMENT", bool is_partitioned=false)
DirichletTemperatureBc(MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x, Vec f)
MoFEMErrorCode iNitialize()
DirichletTemperatureBc(MoFEM::Interface &m_field, const std::string &field_name)
this struct keeps basic methods for moab meshset about material and boundary conditions
Deprecated interface functions.
Data structure to exchange data between mofem and User Loop Methods on entities.
boost::shared_ptr< Field > fieldPtr
boost::shared_ptr< FieldEntity > entPtr
structure for User Loop Methods on finite elements
calculate reactions from vector of internal forces on meshsets
std::string problemName
const ReactionsMap & getReactionsMap() const
Get the Reactions Map.
MoFEM::Interface & mField
ReactionsMap reactionsMap
MoFEMErrorCode calculateReactions(Vec &internal)
calculate reactions from a given vector
std::string fieldName
Reactions(MoFEM::Interface &m_field, string problem_name, string field_name)
std::map< int, VectorDouble > ReactionsMap
const VectorDouble & getReactionsFromSet(const int &id) const
Get the Reactions at specified meshset id.