v0.13.0
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 /* This file is part of MoFEM.
22  * MoFEM is free software: you can redistribute it and/or modify it under
23  * the terms of the GNU Lesser General Public License as published by the
24  * Free Software Foundation, either version 3 of the License, or (at your
25  * option) any later version.
26  *
27  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
28  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
29  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
30  * License for more details.
31  *
32  * You should have received a copy of the GNU Lesser General Public
33  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
34 
35 #ifndef __DIRICHLET_HPP__
36 #define __DIRICHLET_HPP__
37 
38 using namespace boost::numeric;
39 
40 /** \brief Data from Cubit blocksets
41  * \ingroup Dirichlet_bc
42  */
43 struct DataFromBc {
47  Range bc_ents[3];
48 
49  // for rotation
53  double theta;
54 
56  : scaled_values(3), initial_values(3), bc_flags(3), is_rotation(false) {}
57 
58  MoFEMErrorCode getBcData(DisplacementCubitBcData &mydata,
59  const MoFEM::CubitMeshSets *it);
60  MoFEMErrorCode getBcData(std::vector<double> &mydata,
61  const MoFEM::CubitMeshSets *it);
62  MoFEMErrorCode getEntitiesFromBc(MoFEM::Interface &mField,
63  const MoFEM::CubitMeshSets *it);
64 };
65 
66 /** \brief Set Dirichlet boundary conditions on displacements
67  * \ingroup Dirichlet_bc
68  */
70 
72  const std::string fieldName; ///< field name to set Dirichlet BC
73  double dIag; ///< diagonal value set on zeroed column and rows
74 
76  const std::string &field_name, Mat Aij, Vec X, Vec F,
77  string blockset_name = "DISPLACEMENT");
79  const std::string &field_name,
80  string blockset_name = "DISPLACEMENT");
81 
82  std::map<DofIdx, FieldData> mapZeroRows;
83  std::vector<int> dofsIndices;
84  std::vector<double> dofsValues;
85  std::vector<double> dofsXValues;
86  const std::string blocksetName;
87 
88  boost::ptr_vector<MethodForForceScaling> methodsOp;
89  virtual MoFEMErrorCode iNitialize();
90 
91  MoFEMErrorCode preProcess();
92  MoFEMErrorCode operator()() { return 0; }
93  MoFEMErrorCode postProcess();
94  /**
95  * @brief Get the Bc Data From Sets And Blocks object
96  * Use DISPLACEMENT blockset name (default)
97  * with 6 atributes:
98  * 1,2,3 are values of displacements x,y,z
99  * 4,5,6 are flags for x,y,z (0 or 1)
100  * @param bc_data
101  * @return MoFEMErrorCode
102  */
103  MoFEMErrorCode getBcDataFromSetsAndBlocks(std::vector<DataFromBc> &bc_data);
104  /**
105  * @brief Get the Rotation Bc From Block object
106  * Use ROTATION blockset name
107  * with 7 atributes:
108  * 1,2,3 are x,y,z coords of the center of rotation
109  * 4,5,6 are are angular velocities in x,y,z
110  * @param bc_data
111  * @return MoFEMErrorCode
112  */
113  MoFEMErrorCode getRotationBcFromBlock(std::vector<DataFromBc> &bc_data);
114 
115  /**
116  * @brief Calculate displacements from rotation for particular dof
117  * @param dof
118  * @param bc_data
119  * @return MoFEMErrorCode
120  */
121  MoFEMErrorCode calculateRotationForDof(VectorDouble3 &coords,
122  DataFromBc &bc_data);
123  MoFEMErrorCode calculateRotationForDof(EntityHandle ent,
124  DataFromBc &bc_data);
125  MoFEMErrorCode applyScaleBcData(DataFromBc &bc_data);
126 };
127 
132  DataFromBc &data_from_dirichlet_bc)
133  : dirichletBcPtr(dirichlet_bc_ptr), dataFromDirichletBc(data_from_dirichlet_bc) {}
134 
135  MoFEMErrorCode preProcess() { return 0; }
136  MoFEMErrorCode postProcess() { return 0; }
139  auto &mField = dirichletBcPtr->mField;
140  auto &bc_it = dataFromDirichletBc;
141 
142  EntityHandle v = entPtr->getEnt();
143  int coeff = fieldPtr->getNbOfCoeffs();
144  CHKERR dirichletBcPtr->calculateRotationForDof(v, bc_it);
145  for (int i = 0; i != coeff; i++) {
146  if (bc_it.bc_flags[i]) {
147  if (entPtr->getEntType() == MBVERTEX) {
148  entPtr->getEntFieldData()[i] = bc_it.scaled_values[i];
149  } else {
150  entPtr->getEntFieldData()[i] = 0;
151  }
152  }
153  }
154 
156  }
157 };
158 
160  // using BcEntMethodDisp::BcEntMethodDisp;
163  DataFromBc &data_from_dirichlet_bc,
164  string material_positions)
165  : BcEntMethodDisp(dirichlet_bc_ptr, data_from_dirichlet_bc),
166  materialPositions(material_positions) {}
167 
170  EntityHandle ent = entPtr->getEnt();
171  auto &mField = dirichletBcPtr->mField;
172  auto &bc_it = dataFromDirichletBc;
173  EntityHandle v = entPtr->getEnt();
174 
175  const FieldEntity_multiIndex *field_ents;
176  CHKERR mField.get_field_ents(&field_ents);
177  auto &field_ents_by_uid = field_ents->get<Unique_mi_tag>();
178 
179  auto get_coords = [&]() {
180  VectorDouble3 coords(3);
181  if (entPtr->getEntType() == MBVERTEX) {
182  auto eit =
183  field_ents_by_uid.find(FieldEntity::getLocalUniqueIdCalculate(
184  mField.get_field_bit_number(materialPositions), ent));
185  if (eit != field_ents_by_uid.end())
186  noalias(coords) = (*eit)->getEntFieldData();
187  else
188  CHKERR mField.get_moab().get_coords(&ent, 1, &*coords.data().begin());
189  }
190  return coords;
191  };
192 
193  int coeff = fieldPtr->getNbOfCoeffs();
194  auto coords = get_coords();
195 
196  CHKERR dirichletBcPtr->calculateRotationForDof(v, bc_it);
197  for (int i = 0; i != coeff; i++) {
198  if (bc_it.bc_flags[i]) {
199  if (entPtr->getEntType() == MBVERTEX) {
200  entPtr->getEntFieldData()[i] = coords(i) + bc_it.scaled_values[i];
201  } else {
202  entPtr->getEntFieldData()[i] = 0;
203  }
204  }
205  }
206 
208  }
209 };
210 
211 /// \deprecated use DirichletDisplacementBc
213 
214 /** \brief Set Dirichlet boundary conditions on spatial displacements
215  * \ingroup Dirichlet_bc
216  */
218 
220  MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x,
221  Vec f, const std::string material_positions = "MESH_NODE_POSITIONS",
222  const std::string blockset_name = "DISPLACEMENT")
223  : DirichletDisplacementBc(m_field, field_name, aij, x, f, blockset_name),
224  materialPositions(material_positions) {}
225 
227  MoFEM::Interface &m_field, const std::string &field_name,
228  const std::string material_positions = "MESH_NODE_POSITIONS",
229  const std::string blockset_name = "DISPLACEMENT")
230  : DirichletDisplacementBc(m_field, field_name, blockset_name),
231  materialPositions(material_positions) {}
232 
233  std::string materialPositions; ///< name of the field with reference material
234  ///< positions
235  std::vector<std::string> fixFields; ///<
236 
238  MoFEMErrorCode iNitialize();
239 };
240 
241 /// \deprecated use DirichletSpatialPositionsBc
244 
246 
248  const std::string &field_name, Mat aij, Vec x, Vec f)
249  : DirichletDisplacementBc(m_field, field_name, aij, x, f) {}
250 
252  const std::string &field_name)
253  : DirichletDisplacementBc(m_field, field_name) {}
254 
255  MoFEMErrorCode iNitialize();
256 };
257 
258 /// \deprecated use DirichletTemperatureBc
260 
261 /** \brief Fix dofs on entities
262  * \ingroup Dirichlet_bc
263  */
265 
266  Range eNts;
267  std::vector<std::string> fieldNames;
269  const std::string field_name, Mat aij, Vec x,
270  Vec f, Range &ents)
271  : DirichletDisplacementBc(m_field, field_name, aij, x, f), eNts(ents) {
272  fieldNames.push_back(fieldName);
273  }
274 
276  const std::string field_name, Range &ents)
277  : DirichletDisplacementBc(m_field, field_name), eNts(ents) {
278  fieldNames.push_back(fieldName);
279  }
280 
281  MoFEMErrorCode iNitialize();
282  MoFEMErrorCode preProcess();
283  MoFEMErrorCode postProcess();
284 };
285 
286 /** \brief Set Dirichlet boundary conditions on displacements by removing dofs
287  * \ingroup Dirichlet_bc
288  */
290 
291  boost::shared_ptr<vector<DataFromBc>> bcDataPtr;
293  string problemName;
294 
296  const std::string &field_name,
297  const std::string &problem_name,
298  string blockset_name = "DISPLACEMENT",
299  bool is_partitioned = false)
300  : DirichletDisplacementBc(m_field, field_name, blockset_name),
301  problemName(problem_name), isPartitioned(is_partitioned) {}
302 
303  MoFEMErrorCode iNitialize();
304 
305  boost::shared_ptr<EntityMethod> getEntMethodPtr(DataFromBc &data) {
306  return boost::make_shared<BcEntMethodDisp>(this, data);
307  }
308 
309  MoFEMErrorCode preProcess();
310  MoFEMErrorCode operator()() { return 0; }
311  MoFEMErrorCode postProcess() { return 0; }
312 };
313 
314 /** \brief Set Dirichlet boundary conditions on spatial positions by removing dofs
315  * \ingroup Dirichlet_bc
316  */
318 
319  std::string materialPositions;
320 
322  MoFEM::Interface &m_field, const std::string &field_name,
323  const std::string &problem_name,
324  const std::string material_positions = "MESH_NODE_POSITIONS",
325  string blockset_name = "DISPLACEMENT", bool is_partitioned = false)
326  : DirichletDisplacementRemoveDofsBc(m_field, field_name, problem_name, blockset_name,
327  is_partitioned),
328  materialPositions(material_positions) {}
329 
330  boost::shared_ptr<EntityMethod> getEntMethodPtr(DataFromBc &data) {
331  return boost::make_shared<BcEntMethodSpatial>(this, data,
332  materialPositions);
333  }
334 };
335 
336 /// \deprecated use DirichletFixFieldAtEntitiesBc
338 
339 /**
340  * \brief Add boundary conditions form block set having 6 attributes
341  *
342  * First 3 values are magnitudes of dofs e.g. in x,y,z direction and next 3 are
343  flags, respectively.
344  * If flag is false ( = 0), particular dof is not taken into account.
345  Usage in Cubit for displacement:
346  block 1 tri 28 32
347  block 1 name "DISPLACEMENT_1"
348  block 1 attribute count 6
349  block 1 attribute index 1 97 # any value
350  block 1 attribute index 2 0
351  block 1 attribute index 3 0
352  block 1 attribute index 4 0 # flag for x dir
353  block 1 attribute index 5 1 # flag for y dir
354  block 1 attribute index 6 1 # flag for z dir
355  This means that we set zero displacement on y and z dir and on x set
356  direction freely. (value 97 is irrelevant because flag for 1 value is 0
357  (false)) It can be usefull if we want to set boundary conditions directly to
358  triangles e.g, since standard boundary conditions in Cubit allow only using
359  nodeset or surface which might not work with mesh based on facet engine (e.g.
360  STL file)
361  */
363 
365  const std::string &field_name,
366  const std::string &blockset_name, Mat aij,
367  Vec x, Vec f)
368  : DirichletDisplacementBc(m_field, field_name, aij, x, f, blockset_name) {
369  }
370 
372  const std::string &field_name,
373  const std::string &blockset_name)
374  : DirichletDisplacementBc(m_field, field_name, blockset_name) {}
375 };
376 
377 /// \deprecated use DirichletSetFieldFromBlockWithFlags
380 
381 /// \deprecated use DirichletSetFieldFromBlockWithFlags
384 
385 /// \deprecated use DirichletSetFieldFromBlockWithFlags
388 /**
389  * @brief calculate reactions from vector of internal forces on meshsets
390  *
391  * example usage
392  *
393  * \code
394  Vec F_int;
395  DMCreateGlobalVector_MoFEM(dm, &F_int);
396 
397  feRhs->snes_ctx = FEMethod::CTX_SNESSETFUNCTION;
398  feRhs->snes_f = F_int;
399  DMoFEMLoopFiniteElements(dm, "ELASTIC", feRhs);
400 
401  VecAssemblyBegin(F_int);
402  VecAssemblyEnd(F_int);
403  VecGhostUpdateBegin(F_int, INSERT_VALUES, SCATTER_FORWARD);
404  VecGhostUpdateEnd(F_int, INSERT_VALUES, SCATTER_FORWARD);
405 
406  Reactions my_react(m_field, "DM_ELASTIC", "U");
407  my_react.calculateReactions(F_int);
408  int fix_nodes_meshset_id = 1;
409  cout << my_react.getReactionsFromSet(fix_nodes_meshset_id) << endl;
410 
411 * \endcode
412  */
413 struct Reactions {
414 
415  Reactions(MoFEM::Interface &m_field, string problem_name, string field_name)
416  : mField(m_field), problemName(problem_name), fieldName(field_name) {}
417 
418  typedef std::map<int, VectorDouble> ReactionsMap;
420  /**
421  * @brief Get the Reactions Map
422  *
423  * @return const ReactionsMap&
424  */
425  inline const ReactionsMap &getReactionsMap() const { return reactionsMap; }
426  /**
427  * @brief Get the Reactions at specified meshset id
428  *
429  * @param id meshset id (from Cubit)
430  * @return const VectorDouble&
431  */
432  inline const VectorDouble &getReactionsFromSet(const int &id) const {
433  return reactionsMap.at(id);
434  }
435  /**
436  * @brief calculate reactions from a given vector
437  *
438  * @param internal forces vector
439  * @return MoFEMErrorCode
440  */
441  MoFEMErrorCode calculateReactions(Vec &internal);
442 
443 private:
444  std::string problemName;
445  std::string fieldName;
447 };
448 
449 #endif //__DIRICHLET_HPP__
450 
451 /**
452  * \defgroup Dirichlet_bc Dirichlet boundary conditions
453  * \ingroup user_modules
454  **/
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:359
#define DEPRECATED
Definition: definitions.h:30
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
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
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
UBlasVector< int > VectorInt
Definition: Types.hpp:78
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
BcEntMethodSpatial(DirichletDisplacementBc *dirichlet_bc_ptr, DataFromBc &data_from_dirichlet_bc, string material_positions)
MoFEMErrorCode operator()()
function is run for every finite element
Data from Cubit blocksets.
Definition: DirichletBC.hpp:43
VectorInt bc_flags
Definition: DirichletBC.hpp:46
FTensor::Tensor1< double, 3 > t_centr
Definition: DirichletBC.hpp:52
VectorDouble scaled_values
Definition: DirichletBC.hpp:44
FTensor::Tensor1< double, 3 > t_normal
Definition: DirichletBC.hpp:51
VectorDouble initial_values
Definition: DirichletBC.hpp:45
double theta
Definition: DirichletBC.hpp:53
bool is_rotation
Definition: DirichletBC.hpp:50
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
DirichletDisplacementBc(MoFEM::Interface &m_field, const std::string &field_name, Mat Aij, Vec X, Vec F, string blockset_name="DISPLACEMENT")
const std::string blocksetName
Definition: DirichletBC.hpp:86
DirichletDisplacementBc(MoFEM::Interface &m_field, const std::string &field_name, string blockset_name="DISPLACEMENT")
std::map< DofIdx, FieldData > mapZeroRows
Definition: DirichletBC.hpp:82
const std::string fieldName
field name to set Dirichlet BC
Definition: DirichletBC.hpp:72
MoFEMErrorCode operator()()
function is run for every finite element
Definition: DirichletBC.hpp:92
double dIag
diagonal value set on zeroed column and rows
Definition: DirichletBC.hpp:73
std::vector< int > dofsIndices
Definition: DirichletBC.hpp:83
std::vector< double > dofsXValues
Definition: DirichletBC.hpp:85
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: DirichletBC.hpp:88
std::vector< double > dofsValues
Definition: DirichletBC.hpp:84
MoFEM::Interface & mField
Definition: DirichletBC.hpp:71
MoFEMErrorCode calculateRotationForDof(VectorDouble3 &coords, DataFromBc &bc_data)
Calculate displacements from rotation for particular dof.
Set Dirichlet boundary conditions on displacements by removing dofs.
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)
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.
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.
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)
boost::shared_ptr< EntityMethod > getEntMethodPtr(DataFromBc &data)
DirichletTemperatureBc(MoFEM::Interface &m_field, const std::string &field_name, Mat aij, Vec x, Vec f)
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.
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
const VectorDouble & getReactionsFromSet(const int &id) const
Get the Reactions at specified meshset id.
ReactionsMap reactionsMap
std::string fieldName
Reactions(MoFEM::Interface &m_field, string problem_name, string field_name)
std::map< int, VectorDouble > ReactionsMap