v0.9.1
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  // get tangent vector array
283  double *f_tangent_front_mesh_array;
284  CHKERR VecGetArray(smootherData.tangentFrontF,
285  &f_tangent_front_mesh_array);
286  // iterate nodes on tet
287  for (int nn = 0; nn < 4; nn++) {
288 
289  // get indices with Lagrange multiplier at node nn
290  FENumeredDofEntityByNameAndEnt::iterator dit, hi_dit;
291  dit = getFEMethod()
292  ->rowPtr->get<Composite_Name_And_Ent_mi_tag>()
293  .lower_bound(boost::make_tuple(
295  hi_dit = getFEMethod()
296  ->rowPtr->get<Composite_Name_And_Ent_mi_tag>()
297  .upper_bound(boost::make_tuple(
299 
300  // continue if Lagrange are on element
301  if (std::distance(dit, hi_dit) > 0) {
302 
303  FENumeredDofEntityByNameAndEnt::iterator diit, hi_diit;
304 
305  // get mesh node positions at node nn
306  diit = getFEMethod()
307  ->rowPtr->get<Composite_Name_And_Ent_mi_tag>()
308  .lower_bound(boost::make_tuple("MESH_NODE_POSITIONS",
309  getConn()[nn]));
310  hi_diit = getFEMethod()
311  ->rowPtr->get<Composite_Name_And_Ent_mi_tag>()
312  .upper_bound(boost::make_tuple("MESH_NODE_POSITIONS",
313  getConn()[nn]));
314 
315  // iterate over dofs on node nn
316  for (; diit != hi_diit; diit++) {
317  // iterate overt dofs in element column
318  for (int ddd = 0; ddd < nb_col; ddd++) {
319  // check consistency, node has to be at crack front
320  if (rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()] !=
321  diit->get()->getPetscGlobalDofIdx()) {
322  SETERRQ2(
323  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
324  "data inconsistency %d != %d",
325  rowFrontIndices[3 * nn + diit->get()->getDofCoeffIdx()],
326  diit->get()->getPetscGlobalDofIdx());
327  }
328  // dof is not on this partition
329  if (diit->get()->getPetscLocalDofIdx() == -1) {
330  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
331  "data inconsistency");
332  }
333  double g =
334  f_tangent_front_mesh_array[diit->get()
335  ->getPetscLocalDofIdx()] *
336  k(3 * nn + diit->get()->getDofCoeffIdx(), ddd);
337  int lambda_idx = dit->get()->getPetscGlobalDofIdx();
338  CHKERR MatSetValues(getFEMethod()->snes_B, 1, &lambda_idx, 1,
339  &col_indices_ptr[ddd], &g, ADD_VALUES);
340  }
341  }
342  }
343  }
344  CHKERR VecRestoreArray(smootherData.tangentFrontF,
345  &f_tangent_front_mesh_array);
346  }
347 
349  }
350  };
351 };
352 
353 #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:482
ublas::vector< int > rowFrontIndices
Definition: Smoother.hpp:244
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:549
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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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 MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
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:601
MyVolumeFE & feRhs
calculate right hand side for tetrahedral elements
Definition: Smoother.hpp:133
MyVolumeFE & getLoopFeLhs()
get lhs volume element
Definition: Smoother.hpp:136
ublas::vector< boost::shared_ptr< const FEDofEntity >, DofsAllocator > VectorDofs
SNESContext snes_ctx
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
SmootherBlockData & smootherData
Definition: Smoother.hpp:228
data for calculation het 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
boost::shared_ptr< const FENumeredDofEntity_multiIndex > rowPtr
Pointer to finite element rows dofs view.
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