v0.14.0
Smoother.hpp
Go to the documentation of this file.
1 
2 
3 #ifndef __SMOOTHER_HPP__
4 #define __SMOOTHER_HPP__
5 
6 #ifndef WITH_ADOL_C
7 #error "MoFEM need to be compiled with ADOL-C"
8 #endif
9 
10 struct Smoother {
11 
13 
14  bool sTabilised;
17  bool ownVectors;
18 
20  : sTabilised(false), frontF(PETSC_NULL), tangentFrontF(PETSC_NULL),
21  ownVectors(false) {
22  ierr = getOptions();
23  CHKERRABORT(PETSC_COMM_SELF, ierr);
24  }
25 
28  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
29  "Get stabilisation element options", "none");
30  CHKERRG(ierr);
31  PetscBool smoothing_on = sTabilised ? PETSC_TRUE : PETSC_FALSE;
32  CHKERR PetscOptionsBool("-smoothing_stabilise",
33  "all nodes controlled by smoothing element", "",
34  smoothing_on, &smoothing_on, PETSC_NULL);
35  sTabilised = (smoothing_on == PETSC_TRUE) ? true : false;
36  ierr = PetscOptionsEnd();
37  CHKERRG(ierr);
39  }
40 
41  virtual ~SmootherBlockData() {
42  if (ownVectors) {
43  ierr = VecDestroy(&frontF);
44  CHKERRABORT(PETSC_COMM_WORLD, ierr);
45  ierr = VecDestroy(&tangentFrontF);
46  CHKERRABORT(PETSC_COMM_WORLD, ierr);
47  }
48  }
49  };
51 
52  std::map<int, NonlinearElasticElement::BlockData> setOfBlocks;
54 
56 
58 
59  MyVolumeFE(MoFEM::Interface &m_field, SmootherBlockData &smoother_data)
61  smootherData(smoother_data) {}
62 
65 
66  CHKERR VolumeElementForcesAndSourcesCore::preProcess();
67 
68  if (A != PETSC_NULL)
69  snes_B = A;
70 
71  if (F != PETSC_NULL)
72  snes_f = F;
73  switch (snes_ctx) {
74  case CTX_SNESSETFUNCTION: {
75  if (smootherData.frontF) {
76  CHKERR VecZeroEntries(smootherData.frontF);
77  CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
78  SCATTER_FORWARD);
79  CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
80  SCATTER_FORWARD);
81  }
82  } break;
83  default:
84  break;
85  }
86 
88  }
89 
92 
93  switch (snes_ctx) {
94  case CTX_SNESSETFUNCTION: {
95  if (smootherData.frontF) {
96  CHKERR VecAssemblyBegin(smootherData.frontF);
97  CHKERR VecAssemblyEnd(smootherData.frontF);
98  CHKERR VecGhostUpdateBegin(smootherData.frontF, ADD_VALUES,
99  SCATTER_REVERSE);
100  CHKERR VecGhostUpdateEnd(smootherData.frontF, ADD_VALUES,
101  SCATTER_REVERSE);
102  CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
103  SCATTER_FORWARD);
104  CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
105  SCATTER_FORWARD);
106  }
107  break;
108  default:
109  break;
110  }
111  }
112 
113  CHKERR VolumeElementForcesAndSourcesCore::postProcess();
115  }
116  };
117 
118  boost::shared_ptr<MyVolumeFE> feRhsPtr;
119  boost::shared_ptr<MyVolumeFE> feLhsPtr;
120 
121  MyVolumeFE &feRhs; ///< calculate right hand side for tetrahedral elements
122  MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
123  MyVolumeFE &feLhs; //< calculate left hand side for tetrahedral elements
124  MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
125 
127  : feRhsPtr(new MyVolumeFE(m_field, smootherData)),
128  feLhsPtr(new MyVolumeFE(m_field, smootherData)), feRhs(*feRhsPtr),
129  feLhs(*feLhsPtr) {}
130 
133 
134  OpJacobianSmoother(const std::string field_name,
137  int tag, bool jacobian)
139  field_name, data, common_data, tag, jacobian, false, false) {}
140 
143 
144  CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
146 
147  commonData.sTress[gg].resize(3, 3, false);
148  for (int dd1 = 0; dd1 < 3; dd1++) {
149  for (int dd2 = 0; dd2 < 3; dd2++) {
150  dAta.materialAdoublePtr->P(dd1, dd2) >>=
151  (commonData.sTress[gg])(dd1, dd2);
152  }
153  }
154 
156  }
157  };
158 
160 
162 
163  OpRhsSmoother(const std::string field_name,
166  SmootherBlockData &smoother_data)
168  common_data),
169  smootherData(smoother_data) {}
170 
171  ublas::vector<int> frontIndices;
172 
173  MoFEMErrorCode aSemble(int row_side, EntityType row_type,
174  EntitiesFieldData::EntData &row_data) {
176 
177  int nb_dofs = row_data.getIndices().size();
178  int *indices_ptr = &row_data.getIndices()[0];
179 
180  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
181  iNdices.resize(nb_dofs, false);
182  noalias(iNdices) = row_data.getIndices();
183  if (!smootherData.sTabilised) {
184  indices_ptr = &iNdices[0];
185  }
186  frontIndices.resize(nb_dofs, false);
187  noalias(frontIndices) = row_data.getIndices();
188  VectorDofs &dofs = row_data.getFieldDofs();
189  VectorDofs::iterator dit = dofs.begin();
190  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
191  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
193  iNdices[ii] = -1;
194  } else {
195  frontIndices[ii] = -1;
196  }
197  }
198  if (smootherData.frontF) {
200  &nf[0], ADD_VALUES);
201  }
202  }
203 
204  CHKERR VecSetOption(getFEMethod()->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
205  PETSC_TRUE);
206  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
207  ADD_VALUES);
208 
210  }
211  };
212 
215 
218 
220  const std::string vel_field, const std::string field_name,
223  SmootherBlockData &smoother_data,
224  const std::string
225  crack_area_tangent_constrains // = "LAMBDA_CRACK_TANGENT_CONSTRAIN"
226  )
228  data, common_data),
229  smootherData(smoother_data),
230  fieldCrackAreaTangentConstrains(crack_area_tangent_constrains) {}
231 
232  ublas::vector<int> rowFrontIndices;
233 
234  MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
235  EntityType col_type,
236  EntitiesFieldData::EntData &row_data,
237  EntitiesFieldData::EntData &col_data) {
239 
240  int nb_row = row_data.getIndices().size();
241  int nb_col = col_data.getIndices().size();
242  int *row_indices_ptr = &row_data.getIndices()[0];
243  int *col_indices_ptr = &col_data.getIndices()[0];
244 
245  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
246  rowIndices.resize(nb_row, false);
247  noalias(rowIndices) = row_data.getIndices();
248  if (!smootherData.sTabilised) {
249  row_indices_ptr = &rowIndices[0];
250  }
251  rowFrontIndices.resize(nb_row, false);
252  noalias(rowFrontIndices) = row_data.getIndices();
253  VectorDofs &dofs = row_data.getFieldDofs();
254  VectorDofs::iterator dit = dofs.begin();
255  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
256  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
258  rowIndices[ii] = -1;
259  } else {
260  rowFrontIndices[ii] = -1;
261  }
262  }
263  }
264 
265  CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr,
266  nb_col, col_indices_ptr, &k(0, 0), ADD_VALUES);
267 
269 
270  const auto bit_number_for_crack_area_tangent_constrain =
272  const auto bit_number_for_mesh_position =
273  getFEMethod()->getFieldBitNumber("MESH_NODE_POSITIONS");
274 
275  // get tangent vector array
276  double *f_tangent_front_mesh_array;
277  CHKERR VecGetArray(smootherData.tangentFrontF,
278  &f_tangent_front_mesh_array);
279 
280  auto row_dofs = getFEMethod()->getRowDofsPtr();
281 
282  // iterate nodes on tet
283  for (int nn = 0; nn < 4; nn++) {
284 
285  // get indices with Lagrange multiplier at node nn
286  auto dit = row_dofs->get<Unique_mi_tag>().lower_bound(
287  DofEntity::getLoFieldEntityUId(
288  bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
289  auto hi_dit = row_dofs->get<Unique_mi_tag>().upper_bound(
290  DofEntity::getHiFieldEntityUId(
291  bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
292 
293  // continue if Lagrange are on element
294  if (std::distance(dit, hi_dit) > 0) {
295 
296  // get mesh node positions at node nn
297  auto diit = row_dofs->get<Unique_mi_tag>().lower_bound(
298  DofEntity::getLoFieldEntityUId(bit_number_for_mesh_position,
299  getConn()[nn]));
300 
301  auto hi_diit = row_dofs->get<Unique_mi_tag>().upper_bound(
302  DofEntity::getHiFieldEntityUId(bit_number_for_mesh_position,
303  getConn()[nn]));
304 
305  // iterate over dofs on node nn
306  for (; diit != hi_diit; diit++) {
307  // iterate overt dofs in element column
308  for (int ddd = 0; ddd < nb_col; ddd++) {
309  // check consistency, node has to be at crack front
310  if (rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()] !=
311  diit->get()->getPetscGlobalDofIdx()) {
312  SETERRQ2(
313  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
314  "data inconsistency %d != %d",
315  rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()],
316  diit->get()->getPetscGlobalDofIdx());
317  }
318  // dof is not on this partition
319  if (diit->get()->getPetscLocalDofIdx() == -1) {
320  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
321  "data inconsistency");
322  }
323  double g =
324  f_tangent_front_mesh_array[diit->get()
325  ->getPetscLocalDofIdx()] *
326  k(3 * nn + diit->get()->getDofCoeffIdx(), ddd);
327  int lambda_idx = dit->get()->getPetscGlobalDofIdx();
328  CHKERR MatSetValues(getFEMethod()->snes_B, 1, &lambda_idx, 1,
329  &col_indices_ptr[ddd], &g, ADD_VALUES);
330  }
331  }
332  }
333  }
334  CHKERR VecRestoreArray(smootherData.tangentFrontF,
335  &f_tangent_front_mesh_array);
336  }
337 
339  }
340  };
341 };
342 
343 #endif //__SMOOTHER_HPP__
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::k
MatrixDouble k
Definition: NonLinearElasticElement.hpp:571
g
constexpr double g
Definition: shallow_wave.cpp:63
Smoother::SmootherBlockData::frontF
Vec frontF
Definition: Smoother.hpp:15
Smoother::OpLhsSmoother::OpLhsSmoother
OpLhsSmoother(const std::string vel_field, const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, SmootherBlockData &smoother_data, const std::string crack_area_tangent_constrains)
Definition: Smoother.hpp:219
Smoother::OpLhsSmoother
Definition: Smoother.hpp:213
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
NonlinearElasticElement::OpRhsPiolaKirchhoff::dAta
BlockData & dAta
Definition: NonLinearElasticElement.hpp:523
NonlinearElasticElement::BlockData::materialAdoublePtr
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr
Definition: NonLinearElasticElement.hpp:92
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getConn
const EntityHandle * getConn()
get element connectivity
Definition: VolumeElementForcesAndSourcesCore.hpp:156
Smoother::SmootherBlockData::~SmootherBlockData
virtual ~SmootherBlockData()
Definition: Smoother.hpp:41
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::SnesMethod::snes_ctx
SNESContext snes_ctx
Definition: LoopMethods.hpp:118
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
Smoother::MyVolumeFE::smootherData
SmootherBlockData & smootherData
Definition: Smoother.hpp:57
NonlinearElasticElement::MyVolumeFE::A
Mat A
Definition: NonLinearElasticElement.hpp:33
Smoother::OpRhsSmoother::OpRhsSmoother
OpRhsSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, SmootherBlockData &smoother_data)
Definition: Smoother.hpp:163
NonlinearElasticElement::OpRhsPiolaKirchhoff::OpRhsPiolaKirchhoff
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:595
Smoother::feLhsPtr
boost::shared_ptr< MyVolumeFE > feLhsPtr
Definition: Smoother.hpp:119
Smoother::feRhsPtr
boost::shared_ptr< MyVolumeFE > feRhsPtr
Definition: Smoother.hpp:118
Smoother::MyVolumeFE::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: Smoother.hpp:63
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
Smoother::SmootherBlockData::ownVectors
bool ownVectors
Definition: Smoother.hpp:17
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::SnesMethod::CTX_SNESSETFUNCTION
@ CTX_SNESSETFUNCTION
Definition: LoopMethods.hpp:107
NonlinearElasticElement::BlockData::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: HookeElement.hpp:37
Smoother::feLhs
MyVolumeFE & feLhs
Definition: Smoother.hpp:123
Smoother::SmootherBlockData::sTabilised
bool sTabilised
Definition: Smoother.hpp:14
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::OpLhsPiolaKirchhoff_dx
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:723
MoFEM::FEMethod::getRowDofsPtr
auto getRowDofsPtr() const
Definition: LoopMethods.hpp:439
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::OpJacobianPiolaKirchhoffStress
OpJacobianPiolaKirchhoffStress(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian, bool ale, bool field_disp)
Construct operator to calculate Piola-Kirchhoff stress or its derivatives over gradient deformation.
Definition: NonLinearElasticElement.cpp:202
Smoother::MyVolumeFE::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: Smoother.hpp:90
Smoother::OpRhsSmoother::smootherData
SmootherBlockData & smootherData
Definition: Smoother.hpp:161
Smoother::MyVolumeFE
Definition: Smoother.hpp:55
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::commonData
CommonData & commonData
Definition: NonLinearElasticElement.hpp:375
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::rowIndices
ublas::vector< int > rowIndices
Definition: NonLinearElasticElement.hpp:564
Smoother::OpLhsSmoother::fieldCrackAreaTangentConstrains
const std::string fieldCrackAreaTangentConstrains
Definition: Smoother.hpp:217
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
Smoother::OpJacobianSmoother::OpJacobianSmoother
OpJacobianSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, int tag, bool jacobian)
Definition: Smoother.hpp:134
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
Smoother::commonData
NonlinearElasticElement::CommonData commonData
Definition: Smoother.hpp:53
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1000
Smoother::OpLhsSmoother::aSemble
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: Smoother.hpp:234
NonlinearElasticElement::OpRhsPiolaKirchhoff::nf
VectorDouble nf
Definition: NonLinearElasticElement.hpp:532
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
Smoother::OpRhsSmoother::frontIndices
ublas::vector< int > frontIndices
Definition: Smoother.hpp:171
NonlinearElasticElement::MyVolumeFE::F
Vec F
Definition: NonLinearElasticElement.hpp:34
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1042
Smoother::SmootherBlockData::tangentFrontF
Vec tangentFrontF
Definition: Smoother.hpp:16
NonlinearElasticElement::CommonData::sTress
std::vector< MatrixDouble3by3 > sTress
Definition: NonLinearElasticElement.hpp:111
NonlinearElasticElement::MyVolumeFE
definition of volume element
Definition: NonLinearElasticElement.hpp:31
Smoother::OpJacobianSmoother::calculateStress
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition: Smoother.hpp:141
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
Smoother
Definition: Smoother.hpp:10
Smoother::feRhs
MyVolumeFE & feRhs
calculate right hand side for tetrahedral elements
Definition: Smoother.hpp:121
Smoother::SmootherBlockData
Definition: Smoother.hpp:12
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
MoFEM::BasicMethod::getFieldBitNumber
unsigned int getFieldBitNumber(std::string field_name) const
Definition: LoopMethods.hpp:270
Smoother::OpLhsSmoother::smootherData
SmootherBlockData & smootherData
Definition: Smoother.hpp:216
Smoother::smootherData
SmootherBlockData smootherData
Definition: Smoother.hpp:50
Smoother::SmootherBlockData::getOptions
MoFEMErrorCode getOptions()
Definition: Smoother.hpp:26
Smoother::setOfBlocks
std::map< int, NonlinearElasticElement::BlockData > setOfBlocks
Definition: Smoother.hpp:52
Smoother::OpLhsSmoother::rowFrontIndices
ublas::vector< int > rowFrontIndices
Definition: Smoother.hpp:232
Smoother::MyVolumeFE::MyVolumeFE
MyVolumeFE(MoFEM::Interface &m_field, SmootherBlockData &smoother_data)
Definition: Smoother.hpp:59
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
Smoother::Smoother
Smoother(MoFEM::Interface &m_field)
Definition: Smoother.hpp:126
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
NonlinearElasticElement::OpRhsPiolaKirchhoff::iNdices
ublas::vector< int > iNdices
Definition: NonLinearElasticElement.hpp:528
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress::dAta
BlockData & dAta
Definition: NonLinearElasticElement.hpp:373
Smoother::OpRhsSmoother
Definition: Smoother.hpp:159
Smoother::SmootherBlockData::SmootherBlockData
SmootherBlockData()
Definition: Smoother.hpp:19
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
Smoother::OpRhsSmoother::aSemble
MoFEMErrorCode aSemble(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
Definition: Smoother.hpp:173
Smoother::OpJacobianSmoother
Definition: Smoother.hpp:131
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx::dAta
BlockData & dAta
Definition: NonLinearElasticElement.hpp:559
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
NonlinearElasticElement::CommonData
common data used by volume elements
Definition: NonLinearElasticElement.hpp:105