v0.10.0
Smoother.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #ifndef __SMOOTHER_HPP__
16 #define __SMOOTHER_HPP__
17 
18 #ifndef WITH_ADOL_C
19 #error "MoFEM need to be compiled with ADOL-C"
20 #endif
21 
22 struct Smoother {
23 
25 
26  bool sTabilised;
27  Vec frontF;
29  bool ownVectors;
30 
32  : sTabilised(false), frontF(PETSC_NULL), tangentFrontF(PETSC_NULL),
33  ownVectors(false) {
34  ierr = getOptions();
35  CHKERRABORT(PETSC_COMM_SELF, ierr);
36  }
37 
40  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
41  "Get stabilisation element options", "none");
42  CHKERRG(ierr);
43  PetscBool smoothing_on = sTabilised ? PETSC_TRUE : PETSC_FALSE;
44  CHKERR PetscOptionsBool("-smoothing_stabilise",
45  "all nodes controlled by smoothing element", "",
46  smoothing_on, &smoothing_on, PETSC_NULL);
47  sTabilised = (smoothing_on == PETSC_TRUE) ? true : false;
48  ierr = PetscOptionsEnd();
49  CHKERRG(ierr);
51  }
52 
53  virtual ~SmootherBlockData() {
54  if (ownVectors) {
55  ierr = VecDestroy(&frontF);
56  CHKERRABORT(PETSC_COMM_WORLD, ierr);
57  ierr = VecDestroy(&tangentFrontF);
58  CHKERRABORT(PETSC_COMM_WORLD, ierr);
59  }
60  }
61  };
63 
64  std::map<int, NonlinearElasticElement::BlockData> setOfBlocks;
66 
68 
70 
71  MyVolumeFE(MoFEM::Interface &m_field, SmootherBlockData &smoother_data)
73  smootherData(smoother_data) {}
74 
77 
78  CHKERR VolumeElementForcesAndSourcesCore::preProcess();
79 
80  if (A != PETSC_NULL)
81  snes_B = A;
82 
83  if (F != PETSC_NULL)
84  snes_f = F;
85  switch (snes_ctx) {
86  case CTX_SNESSETFUNCTION: {
87  if (smootherData.frontF) {
88  CHKERR VecZeroEntries(smootherData.frontF);
89  CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
90  SCATTER_FORWARD);
91  CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
92  SCATTER_FORWARD);
93  }
94  } break;
95  default:
96  break;
97  }
98 
100  }
101 
104 
105  switch (snes_ctx) {
106  case CTX_SNESSETFUNCTION: {
107  if (smootherData.frontF) {
108  CHKERR VecAssemblyBegin(smootherData.frontF);
109  CHKERR VecAssemblyEnd(smootherData.frontF);
110  CHKERR VecGhostUpdateBegin(smootherData.frontF, ADD_VALUES,
111  SCATTER_REVERSE);
112  CHKERR VecGhostUpdateEnd(smootherData.frontF, ADD_VALUES,
113  SCATTER_REVERSE);
114  CHKERR VecGhostUpdateBegin(smootherData.frontF, INSERT_VALUES,
115  SCATTER_FORWARD);
116  CHKERR VecGhostUpdateEnd(smootherData.frontF, INSERT_VALUES,
117  SCATTER_FORWARD);
118  }
119  break;
120  default:
121  break;
122  }
123  }
124 
125  CHKERR VolumeElementForcesAndSourcesCore::postProcess();
127  }
128  };
129 
130  boost::shared_ptr<MyVolumeFE> feRhsPtr;
131  boost::shared_ptr<MyVolumeFE> feLhsPtr;
132 
133  MyVolumeFE &feRhs; ///< calculate right hand side for tetrahedral elements
134  MyVolumeFE &getLoopFeRhs() { return feRhs; } ///< get rhs volume element
135  MyVolumeFE &feLhs; //< calculate left hand side for tetrahedral elements
136  MyVolumeFE &getLoopFeLhs() { return feLhs; } ///< get lhs volume element
137 
139  : feRhsPtr(new MyVolumeFE(m_field, smootherData)),
140  feLhsPtr(new MyVolumeFE(m_field, smootherData)), feRhs(*feRhsPtr),
141  feLhs(*feLhsPtr) {}
142 
145 
146  OpJacobianSmoother(const std::string field_name,
149  int tag, bool jacobian)
151  field_name, data, common_data, tag, jacobian, false, false) {}
152 
155 
156  CHKERR dAta.materialAdoublePtr->calculateP_PiolaKirchhoffI(
158 
159  commonData.sTress[gg].resize(3, 3, false);
160  for (int dd1 = 0; dd1 < 3; dd1++) {
161  for (int dd2 = 0; dd2 < 3; dd2++) {
162  dAta.materialAdoublePtr->P(dd1, dd2) >>=
163  (commonData.sTress[gg])(dd1, dd2);
164  }
165  }
166 
168  }
169  };
170 
172 
174 
175  OpRhsSmoother(const std::string field_name,
178  SmootherBlockData &smoother_data)
179  : NonlinearElasticElement::OpRhsPiolaKirchhoff(field_name, data,
180  common_data),
181  smootherData(smoother_data) {}
182 
183  ublas::vector<int> frontIndices;
184 
185  MoFEMErrorCode aSemble(int row_side, EntityType row_type,
188 
189  int nb_dofs = row_data.getIndices().size();
190  int *indices_ptr = &row_data.getIndices()[0];
191 
192  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
193  iNdices.resize(nb_dofs, false);
194  noalias(iNdices) = row_data.getIndices();
195  if (!smootherData.sTabilised) {
196  indices_ptr = &iNdices[0];
197  }
198  frontIndices.resize(nb_dofs, false);
199  noalias(frontIndices) = row_data.getIndices();
200  VectorDofs &dofs = row_data.getFieldDofs();
201  VectorDofs::iterator dit = dofs.begin();
202  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
203  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
205  iNdices[ii] = -1;
206  } else {
207  frontIndices[ii] = -1;
208  }
209  }
210  if (smootherData.frontF) {
212  &nf[0], ADD_VALUES);
213  }
214  }
215 
216  CHKERR VecSetOption(getFEMethod()->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
217  PETSC_TRUE);
218  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nf[0],
219  ADD_VALUES);
220 
222  }
223  };
224 
227 
230 
232  const std::string vel_field, const std::string field_name,
235  SmootherBlockData &smoother_data,
236  const std::string
237  crack_area_tangent_constrains // = "LAMBDA_CRACK_TANGENT_CONSTRAIN"
238  )
239  : NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(vel_field, field_name,
240  data, common_data),
241  smootherData(smoother_data),
242  fieldCrackAreaTangentConstrains(crack_area_tangent_constrains) {}
243 
244  ublas::vector<int> rowFrontIndices;
245 
246  MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type,
247  EntityType col_type,
251 
252  int nb_row = row_data.getIndices().size();
253  int nb_col = col_data.getIndices().size();
254  int *row_indices_ptr = &row_data.getIndices()[0];
255  int *col_indices_ptr = &col_data.getIndices()[0];
256 
257  if (!dAta.forcesOnlyOnEntitiesRow.empty()) {
258  rowIndices.resize(nb_row, false);
259  noalias(rowIndices) = row_data.getIndices();
260  if (!smootherData.sTabilised) {
261  row_indices_ptr = &rowIndices[0];
262  }
263  rowFrontIndices.resize(nb_row, false);
264  noalias(rowFrontIndices) = row_data.getIndices();
265  VectorDofs &dofs = row_data.getFieldDofs();
266  VectorDofs::iterator dit = dofs.begin();
267  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
268  if (dAta.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) !=
270  rowIndices[ii] = -1;
271  } else {
272  rowFrontIndices[ii] = -1;
273  }
274  }
275  }
276 
277  CHKERR MatSetValues(getFEMethod()->snes_B, nb_row, row_indices_ptr,
278  nb_col, col_indices_ptr, &k(0, 0), ADD_VALUES);
279 
281 
282  const auto bit_number_for_crack_area_tangent_constrain =
284  const auto bit_number_for_mesh_position =
285  getFEMethod()->getFieldBitNumber("MESH_NODE_POSITIONS");
286 
287  // get tangent vector array
288  double *f_tangent_front_mesh_array;
289  CHKERR VecGetArray(smootherData.tangentFrontF,
290  &f_tangent_front_mesh_array);
291 
292  auto row_dofs = getFEMethod()->getRowDofsPtr();
293 
294  // iterate nodes on tet
295  for (int nn = 0; nn < 4; nn++) {
296 
297  // get indices with Lagrange multiplier at node nn
298  auto dit = row_dofs->get<Unique_mi_tag>().lower_bound(
299  DofEntity::getLoFieldEntityUId(
300  bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
301  auto hi_dit = row_dofs->get<Unique_mi_tag>().upper_bound(
302  DofEntity::getHiFieldEntityUId(
303  bit_number_for_crack_area_tangent_constrain, getConn()[nn]));
304 
305  // continue if Lagrange are on element
306  if (std::distance(dit, hi_dit) > 0) {
307 
308  // get mesh node positions at node nn
309  auto diit = row_dofs->get<Unique_mi_tag>().lower_bound(
310  DofEntity::getLoFieldEntityUId(bit_number_for_mesh_position,
311  getConn()[nn]));
312 
313  auto hi_diit = row_dofs->get<Unique_mi_tag>().upper_bound(
314  DofEntity::getHiFieldEntityUId(bit_number_for_mesh_position,
315  getConn()[nn]));
316 
317  // iterate over dofs on node nn
318  for (; diit != hi_diit; diit++) {
319  // iterate overt dofs in element column
320  for (int ddd = 0; ddd < nb_col; ddd++) {
321  // check consistency, node has to be at crack front
322  if (rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()] !=
323  diit->get()->getPetscGlobalDofIdx()) {
324  SETERRQ2(
325  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
326  "data inconsistency %d != %d",
327  rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()],
328  diit->get()->getPetscGlobalDofIdx());
329  }
330  // dof is not on this partition
331  if (diit->get()->getPetscLocalDofIdx() == -1) {
332  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
333  "data inconsistency");
334  }
335  double g =
336  f_tangent_front_mesh_array[diit->get()
337  ->getPetscLocalDofIdx()] *
338  k(3 * nn + diit->get()->getDofCoeffIdx(), ddd);
339  int lambda_idx = dit->get()->getPetscGlobalDofIdx();
340  CHKERR MatSetValues(getFEMethod()->snes_B, 1, &lambda_idx, 1,
341  &col_indices_ptr[ddd], &g, ADD_VALUES);
342  }
343  }
344  }
345  }
346  CHKERR VecRestoreArray(smootherData.tangentFrontF,
347  &f_tangent_front_mesh_array);
348  }
349 
351  }
352  };
353 };
354 
355 #endif //__SMOOTHER_HPP__
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.
Deprecated interface functions.
MoFEMErrorCode aSemble(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: Smoother.hpp:185
boost::shared_ptr< FunctionsToCalculatePiolaKirchhoffI< adouble > > materialAdoublePtr
MyVolumeFE(MoFEM::Interface &m_field, SmootherBlockData &smoother_data)
Definition: Smoother.hpp:71
std::vector< MatrixDouble3by3 > sTress
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: Smoother.hpp:102
OpRhsSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, SmootherBlockData &smoother_data)
Definition: Smoother.hpp:175
MyVolumeFE & feLhs
Definition: Smoother.hpp:135
MoFEMErrorCode getOptions()
Definition: Smoother.hpp:38
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
ublas::vector< int > rowFrontIndices
Definition: Smoother.hpp:244
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
MoFEMErrorCode aSemble(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: Smoother.hpp:246
unsigned int getFieldBitNumber(std::string field_name) const
boost::shared_ptr< MyVolumeFE > feLhsPtr
Definition: Smoother.hpp:131
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: Smoother.hpp:75
Vec & snes_f
residual
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
SmootherBlockData smootherData
Definition: Smoother.hpp:62
ublas::vector< int > frontIndices
Definition: Smoother.hpp:183
Mat & snes_B
preconditioner of jacobian matrix
OpLhsPiolaKirchhoff_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
common data used by volume elements
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
DataForcesAndSourcesCore::EntData EntData
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
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:231
NonlinearElasticElement::CommonData commonData
Definition: Smoother.hpp:65
MoFEMErrorCode calculateStress(const int gg)
Calculate Paola-Kirchhoff I stress.
Definition: Smoother.hpp:153
MyVolumeFE & getLoopFeRhs()
get rhs volume element
Definition: Smoother.hpp:134
Smoother(MoFEM::Interface &m_field)
Definition: Smoother.hpp:138
boost::shared_ptr< MyVolumeFE > feRhsPtr
Definition: Smoother.hpp:130
std::map< int, NonlinearElasticElement::BlockData > setOfBlocks
Definition: Smoother.hpp:64
#define CHKERR
Inline error check.
Definition: definitions.h:604
MyVolumeFE & feRhs
calculate right hand side for tetrahedral elements
Definition: Smoother.hpp:133
MyVolumeFE & getLoopFeLhs()
get lhs volume element
Definition: Smoother.hpp:136
SNESContext snes_ctx
auto getRowDofsPtr() const
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
SmootherBlockData & smootherData
Definition: Smoother.hpp:228
data for calculation heat conductivity and heat capacity elements
OpJacobianSmoother(const std::string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, int tag, bool jacobian)
Definition: Smoother.hpp:146
const std::string fieldCrackAreaTangentConstrains
Definition: Smoother.hpp:229
SmootherBlockData & smootherData
Definition: Smoother.hpp:173
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
SmootherBlockData & smootherData
Definition: Smoother.hpp:69