v0.11.1
mesh_smoothing.cpp
Go to the documentation of this file.
1 /** \file mesh_smoothing.cpp
2  * \brief Improve mesh quality using Volume-length quality measure with barrier
3  * \example mesh_smoothing.cpp
4  *
5  */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #include <MoFEM.hpp>
22 using namespace MoFEM;
23 
24 #include <BasicFiniteElements.hpp>
25 
26 #include <Smoother.hpp>
28 #include <VolumeLengthQuality.hpp>
29 
30 using namespace MoFEM;
31 
32 static char help[] = "mesh smoothing\n\n";
33 
34 PetscBool flg_myfile = PETSC_TRUE;
35 char mesh_file_name[255];
38 PetscBool output_vtk = PETSC_TRUE;
39 double tol = 0.1;
40 double gamma_factor = 0.8;
41 
42 int main(int argc, char *argv[]) {
43 
44  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
45 
46  try {
47 
48  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Mesh cut options", "none");
49  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
50  mesh_file_name, 255, &flg_myfile);
51  CHKERR PetscOptionsInt("-edges_block_set", "edges side set", "",
52  edges_block_set, &edges_block_set, PETSC_NULL);
53  CHKERR PetscOptionsInt("-vertex_block_set", "vertex side set", "",
54  vertex_block_set, &vertex_block_set, PETSC_NULL);
55  CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
56  output_vtk, &output_vtk, PETSC_NULL);
57  CHKERR PetscOptionsScalar("-quality_reduction_tol", "",
58  "Tolerance of quality reduction", tol, &tol,
59  PETSC_NULL);
60  CHKERR PetscOptionsScalar("-gamma_factor", "",
61  "Gamma factor", gamma_factor, &gamma_factor,
62  PETSC_NULL);
63  ierr = PetscOptionsEnd();
64  CHKERRG(ierr);
65 
66  // Create MoAB database
67  moab::Core moab_core; // create database
68  moab::Interface &moab = moab_core; // create interface to database
69  // Create MoFEM database and link it to MoAB
70  MoFEM::Core mofem_core(moab); // create database
71  MoFEM::Interface &m_field = mofem_core; // create interface to database
72  // Register DM Manager
73  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
74 
75  // Get simple interface is simplified version enabling quick and
76  // easy construction of problem.
77  Simple *simple_interface;
78  // Query interface and get pointer to Simple interface
79  CHKERR m_field.getInterface(simple_interface);
80 
81  // Build problem with simple interface
82  {
83  // Get options for simple interface from command line
84  CHKERR simple_interface->getOptions();
85  // Load mesh file to database
86  CHKERR simple_interface->loadFile();
87 
88  // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
89  // used to construct base functions.
90  CHKERR simple_interface->addDomainField("MESH_NODE_POSITIONS", H1,
92  // Add Lagrange multiplier field on body boundary
93  CHKERR simple_interface->addBoundaryField("LAMBDA_SURFACE", H1,
95 
96  // Set fields order domain and boundary fields.
97  CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS",
98  1); // to approximate function
99  CHKERR simple_interface->setFieldOrder("LAMBDA_SURFACE",
100  1); // to Lagrange multipliers
101 
102  simple_interface->getDomainFEName() = "SMOOTHING";
103  simple_interface->getBoundaryFEName() = "SURFACE_SLIDING";
104 
105  // Other fields and finite elements added to database directly
106  {
107  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(
109  // Declare approximation fields
110  CHKERR m_field.add_field("LAMBDA_EDGE", H1, AINSWORTH_LOBATTO_BASE, 2,
111  MB_TAG_SPARSE, MF_ZERO);
112  Range edges;
114  ->getEntitiesByDimension(edges_block_set, BLOCKSET, 1, edges,
115  true);
116  CHKERR m_field.add_ents_to_field_by_type(edges, MBEDGE,
117  "LAMBDA_EDGE");
119  ->synchroniseFieldEntities("LAMBDA_EDGE");
120  CHKERR m_field.set_field_order(0, MBVERTEX, "LAMBDA_EDGE", 1);
121 
122  CHKERR m_field.add_finite_element("EDGE_SLIDING");
124  "EDGE_SLIDING", "MESH_NODE_POSITIONS");
126  "EDGE_SLIDING", "MESH_NODE_POSITIONS");
128  "EDGE_SLIDING", "MESH_NODE_POSITIONS");
129  CHKERR m_field.modify_finite_element_add_field_row("EDGE_SLIDING",
130  "LAMBDA_EDGE");
131  CHKERR m_field.modify_finite_element_add_field_col("EDGE_SLIDING",
132  "LAMBDA_EDGE");
133  CHKERR m_field.modify_finite_element_add_field_data("EDGE_SLIDING",
134  "LAMBDA_EDGE");
135  CHKERR m_field.add_ents_to_finite_element_by_type(edges, MBEDGE,
136  "EDGE_SLIDING");
137  simple_interface->getOtherFiniteElements().push_back("EDGE_SLIDING");
138  }
139  }
140 
141  CHKERR simple_interface->defineFiniteElements();
142  CHKERR simple_interface->defineProblem();
143  CHKERR simple_interface->buildFields();
144  // Remove vertices form LAMBDA_SURFACE which are on the edges
145  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(
147  Range edges;
148  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
149  edges_block_set, BLOCKSET, 1, edges, true);
150  Range verts;
151  CHKERR m_field.get_moab().get_connectivity(edges, verts, true);
152  CHKERR m_field.remove_ents_from_field("LAMBDA_SURFACE", verts);
153  }
154  CHKERR simple_interface->buildFiniteElements();
155  CHKERR simple_interface->buildProblem();
156  }
157 
158  struct ElementsAndOperators {
159 
160  MoFEM::Interface &mField;
161  Vec minQualityVec;
162  double *minQualityPtr;
163 
164  ElementsAndOperators(MoFEM::Interface &m_field) : mField(m_field) {
165  ierr = VecCreateMPI(PETSC_COMM_WORLD, 1, m_field.get_comm_size(),
166  &minQualityVec);
167  CHKERRABORT(PETSC_COMM_WORLD, ierr);
168  ierr = VecGetArray(minQualityVec, &minQualityPtr);
169  CHKERRABORT(PETSC_COMM_WORLD, ierr);
170  }
171 
172  virtual ~ElementsAndOperators() {
173  ierr = VecRestoreArray(minQualityVec, &minQualityPtr);
174  CHKERRABORT(PETSC_COMM_WORLD, ierr);
175  ierr = VecDestroy(&minQualityVec);
176  CHKERRABORT(PETSC_COMM_WORLD, ierr);
177  }
178 
179  double getMinQuality() const { return *minQualityPtr; }
180 
181  enum Tags {
182  SMOOTHING_TAG = 1,
183  SURFACE_CONSTRAINS_TAG,
184  EDGE_CONSTRAINS_TAG
185  };
186 
187  boost::shared_ptr<Smoother> smootherFe;
188  boost::shared_ptr<Smoother::MyVolumeFE>
189  feSmootherRhs; ///< Integrate smoothing operators
190  boost::shared_ptr<Smoother::MyVolumeFE>
191  feSmootherLhs; ///< Integrate smoothing operators
192  boost::shared_ptr<VolumeLengthQuality<double> > volumeLengthDouble;
193  boost::shared_ptr<VolumeLengthQuality<adouble> > volumeLengthAdouble;
194 
195  boost::shared_ptr<SurfaceSlidingConstrains> surfaceConstrain;
196  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
197  skinOrientation;
198  boost::shared_ptr<EdgeSlidingConstrains> edgeConstrain;
199 
200  boost::shared_ptr<DirichletFixFieldAtEntitiesBc> fixMaterialEnts;
201 
202  boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore> minQualityFe;
203 
204  struct MinQualityOp : MoFEM::ForcesAndSourcesCore::UserDataOperator {
205  double *minQualityPtr;
206  MinQualityOp(double *min_quality_ptr)
208  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
209  minQualityPtr(min_quality_ptr) {}
210  MoFEMErrorCode doWork(int side, EntityType type,
213  if (type != MBVERTEX)
215  double q = Tools::volumeLengthQuality(&*data.getFieldData().begin());
216  *minQualityPtr = fmin(*minQualityPtr, q);
218  }
219  };
220 
221  MoFEMErrorCode createSmoothingFE() {
223  smootherFe = boost::shared_ptr<Smoother>(new Smoother(mField));
224  volumeLengthAdouble = boost::shared_ptr<VolumeLengthQuality<adouble> >(
226  volumeLengthDouble = boost::shared_ptr<VolumeLengthQuality<double> >(
228 
229  Range tets;
230  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, tets);
231  smootherFe->setOfBlocks[0].tEts.merge(tets);
232 
233  smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
234  smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
235 
236  // set element data
237  smootherFe->commonData.spatialPositions = "MESH_NODE_POSITIONS";
238  smootherFe->commonData.meshPositions = "NONE";
239 
240  smootherFe->feRhs.meshPositionsFieldName = "NONE";
241  smootherFe->feLhs.meshPositionsFieldName = "NONE";
242  smootherFe->feRhs.addToRule = 0;
243  smootherFe->feLhs.addToRule = 0;
244 
245  feSmootherRhs = smootherFe->feRhsPtr;
246  feSmootherLhs = smootherFe->feLhsPtr;
247 
248  // Smoother right hand side
249  feSmootherRhs->getOpPtrVector().push_back(
251  "MESH_NODE_POSITIONS", smootherFe->commonData));
252  feSmootherRhs->getOpPtrVector().push_back(
254  "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
255  smootherFe->commonData, SMOOTHING_TAG, false));
256  feSmootherRhs->getOpPtrVector().push_back(new Smoother::OpRhsSmoother(
257  "MESH_NODE_POSITIONS", smootherFe->setOfBlocks[0],
258  smootherFe->commonData, smootherFe->smootherData));
259 
260  // Smoother left hand side
261  feSmootherLhs->getOpPtrVector().push_back(
263  "MESH_NODE_POSITIONS", smootherFe->commonData));
264  feSmootherLhs->getOpPtrVector().push_back(
266  "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
267  smootherFe->commonData, SMOOTHING_TAG, true));
268  feSmootherLhs->getOpPtrVector().push_back(new Smoother::OpLhsSmoother(
269  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
270  smootherFe->setOfBlocks.at(0), smootherFe->commonData,
271  smootherFe->smootherData, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
272 
273  minQualityFe =
274  boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore>(
276  minQualityFe->getOpPtrVector().push_back(
277  new MinQualityOp(minQualityPtr));
278 
279  Range fixed_vertex;
282  CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
283  vertex_block_set, BLOCKSET, 0, fixed_vertex, true);
284  }
285  fixMaterialEnts = boost::shared_ptr<DirichletFixFieldAtEntitiesBc>(
286  new DirichletFixFieldAtEntitiesBc(mField, "MESH_NODE_POSITIONS",
287  fixed_vertex));
288  fixMaterialEnts->fieldNames.push_back("LAMBDA_SURFACE");
289  fixMaterialEnts->fieldNames.push_back("LAMBDA_EDGE");
290 
292  }
293 
294  MoFEMErrorCode createConstrians() {
296  skinOrientation = boost::shared_ptr<
299  surfaceConstrain = boost::shared_ptr<SurfaceSlidingConstrains>(
300  skinOrientation,
301  new SurfaceSlidingConstrains(mField, *skinOrientation));
302  surfaceConstrain->setOperators(SURFACE_CONSTRAINS_TAG, "LAMBDA_SURFACE",
303  "MESH_NODE_POSITIONS");
304 
307  Range edges;
309  ->getEntitiesByDimension(edges_block_set, BLOCKSET, 1, edges,
310  true);
311 
312  Range tets;
313  CHKERR mField.get_moab().get_entities_by_type(0, MBTET, tets);
314  Skinner skin(&mField.get_moab());
315  Range skin_faces; // skin faces from 3d ents
316  CHKERR skin.find_skin(0, tets, false, skin_faces);
317 
318  edgeConstrain = boost::shared_ptr<EdgeSlidingConstrains>(
319  new EdgeSlidingConstrains(mField));
320  CHKERR edgeConstrain->setOperators(EDGE_CONSTRAINS_TAG, edges,
321  skin_faces, "LAMBDA_EDGE",
322  "MESH_NODE_POSITIONS");
323 
324  // CHKERR EdgeSlidingConstrains::CalculateEdgeBase::saveEdges(
325  // mField.get_moab(), "out_edges.vtk", edges);
326 
327  }
329  }
330 
331  MoFEMErrorCode addFEtoDM(DM dm) {
333  boost::shared_ptr<ForcesAndSourcesCore> null;
334 
335  CHKERR DMMoFEMSNESSetFunction(dm, DM_NO_ELEMENT, null, fixMaterialEnts,
336  null);
337  CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", feSmootherRhs, null,
338  null);
339  CHKERR DMMoFEMSNESSetFunction(dm, "SURFACE_SLIDING",
340  surfaceConstrain->feRhsPtr, null, null);
341  CHKERR DMMoFEMSNESSetFunction(dm, "EDGE_SLIDING",
342  edgeConstrain->feRhsPtr, null, null);
344  fixMaterialEnts);
345 
346  CHKERR DMMoFEMSNESSetJacobian(dm, DM_NO_ELEMENT, null, fixMaterialEnts,
347  null);
348  CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", feSmootherLhs, null,
349  null);
350  CHKERR DMMoFEMSNESSetJacobian(dm, "SURFACE_SLIDING",
351  surfaceConstrain->feLhsPtr, null, null);
352  CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE_SLIDING",
353  edgeConstrain->feLhsPtr, null, null);
355  fixMaterialEnts);
356 
357  // MoFEM::SnesCtx *snes_ctx;
358  // DMMoFEMGetSnesCtx(dm,&snes_ctx);
359  // snes_ctx->vErify = true;
360 
362  }
363 
364  MoFEMErrorCode calcuteMinQuality(DM dm) {
366  *minQualityPtr = 1;
367  CHKERR DMoFEMLoopFiniteElements(dm, "SMOOTHING", minQualityFe.get());
368  CHKERR VecMin(minQualityVec, PETSC_NULL, minQualityPtr);
370  }
371  };
372 
373  ElementsAndOperators elements_and_operators(m_field);
374  CHKERR elements_and_operators.createSmoothingFE();
375  CHKERR elements_and_operators.createConstrians();
376 
377  DM dm;
378  CHKERR simple_interface->getDM(&dm);
379  CHKERR elements_and_operators.addFEtoDM(dm);
380 
381  struct Solve {
382 
383  MoFEMErrorCode operator()(DM dm) const {
385 
386  // Create the right hand side vector and vector of unknowns
387  Vec F, D;
388  CHKERR DMCreateGlobalVector(dm, &F);
389  // Create unknown vector by creating duplicate copy of F vector. only
390  // structure is duplicated no values.
391  CHKERR VecDuplicate(F, &D);
392 
393  CHKERR zeroLambdaFields(dm);
394  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
395 
396  // Create solver and link it to DM
397  SNES solver;
398  CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
399  CHKERR SNESSetFromOptions(solver);
400  CHKERR SNESSetDM(solver, dm);
401  // Set-up solver, is type of solver and pre-conditioners
402  CHKERR SNESSetUp(solver);
403  // At solution process, KSP solver using DM creates matrices, Calculate
404  // values of the left hand side and the right hand side vector. then
405  // solves system of equations. Results are stored in vector D.
406  CHKERR SNESSolve(solver, F, D);
407 
408  // Scatter solution on the mesh. Stores unknown vector on field on the
409  // mesh.
410  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
411  // Clean data. Solver and vector are not needed any more.
412  CHKERR SNESDestroy(&solver);
413  CHKERR VecDestroy(&D);
414  CHKERR VecDestroy(&F);
415 
417  }
418 
419  MoFEMErrorCode setCoordsFromField(DM dm) const {
421  MoFEM::Interface *m_field_ptr;
422  CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
423  for (_IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(*m_field_ptr,
424  "MESH_NODE_POSITIONS", it)) {
425  if (it->get()->getEntType() != MBVERTEX)
426  continue;
427  VectorDouble3 coords(3);
428  for(int dd = 0;dd!=3;++dd)
429  coords[dd] = it->get()->getEntFieldData()[dd];
430  EntityHandle ent = it->get()->getEnt();
431  CHKERR m_field_ptr->get_moab().set_coords(&ent, 1, &*coords.begin());
432  }
434  }
435 
436  MoFEMErrorCode setFieldFromCoords(DM dm) const {
438  MoFEM::Interface *m_field_ptr;
439  CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
440  for (_IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(*m_field_ptr,
441  "MESH_NODE_POSITIONS", it)) {
442  if (it->get()->getEntType() != MBVERTEX)
443  continue;
444  EntityHandle ent = it->get()->getEnt();
445  VectorDouble3 coords(3);
446  CHKERR m_field_ptr->get_moab().get_coords(&ent, 1, &*coords.begin());
447  for(int dd = 0;dd!=3;++dd)
448  it->get()->getEntFieldData()[dd] = coords[dd];
449  }
451  }
452 
453  private:
454  MoFEMErrorCode zeroLambdaFields(DM dm) const {
456  MoFEM::Interface *m_field_ptr;
457  CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
458  CHKERR m_field_ptr->getInterface<FieldBlas>()->setField(
459  0, MBVERTEX, "LAMBDA_SURFACE");
461  }
462 
463  };
464 
465  Solve solve;
466  CHKERR solve.setFieldFromCoords(dm);
467 
468  CHKERR elements_and_operators.calcuteMinQuality(dm);
469  double min_quality = elements_and_operators.getMinQuality();
470  PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f\n", min_quality);
471 
472  double gamma = min_quality > 0 ? gamma_factor * min_quality
473  : min_quality / gamma_factor;
474  elements_and_operators.volumeLengthDouble->gAmma = gamma;
475  elements_and_operators.volumeLengthAdouble->gAmma = gamma;
476 
477  double min_quality_p, eps;
478  do {
479 
480  min_quality_p = min_quality;
481 
482  CHKERR solve(dm);
483 
484  CHKERR solve.setCoordsFromField(dm);
485  CHKERR elements_and_operators.calcuteMinQuality(dm);
486  min_quality = elements_and_operators.getMinQuality();
487 
488  eps = (min_quality - min_quality_p) / min_quality;
489  PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f eps = %4.3f\n",
490  min_quality, eps);
491 
492  double gamma = min_quality > 0 ? gamma_factor * min_quality
493  : min_quality / gamma_factor;
494  elements_and_operators.volumeLengthDouble->gAmma = gamma;
495  elements_and_operators.volumeLengthAdouble->gAmma = gamma;
496 
497  } while (eps > tol);
498 
499  // if (m_field.getInterface<MeshsetsManager>()->checkMeshset(edges_block_set,
500  // BLOCKSET)) {
501  // Range edges;
502  // CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
503  // edges_block_set, BLOCKSET, 1, edges, true);
504 
505  // Range tets;
506  // CHKERR moab.get_entities_by_type(0,MBTET,tets);
507  // Skinner skin(&moab);
508  // Range skin_faces; // skin faces from 3d ents
509  // CHKERR skin.find_skin(0, tets, false, skin_faces);
510 
511  // CHKERR EdgeSlidingConstrains::CalculateEdgeBase::setTags(moab, edges,
512  // skin_faces);
513  // CHKERR EdgeSlidingConstrains::CalculateEdgeBase::saveEdges(
514  // moab, "out_edges.vtk", edges);
515  // }
516 
517  if (output_vtk)
518  CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByType(
519  BitRefLevel().set(0), BitRefLevel().set(), MBTET, "out.vtk", "VTK",
520  "");
521  }
522  CATCH_ERRORS;
523 
525 }
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
Implementing surface sliding constrains.
Implementation of Volume-Length-Quality measure with barrier.
@ BARRIER_AND_QUALITY
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
@ MF_ZERO
Definition: definitions.h:189
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:154
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
@ BLOCKSET
Definition: definitions.h:217
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:637
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:678
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMMoFEM.cpp:336
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:457
#define _IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
loop over all dofs from a moFEM field and particular field
Definition: Interface.hpp:1785
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode remove_ents_from_field(const std::string name, const EntityHandle meshset, const EntityType type, int verb=DEFAULT_VERBOSITY)=0
remove entities from field
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
int main(int argc, char *argv[])
double tol
PetscBool output_vtk
static char help[]
PetscBool flg_myfile
double gamma_factor
int vertex_block_set
int edges_block_set
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1140
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1953
const double D
diffusivity
Fix dofs on entities.
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorDouble & getFieldData() const
get dofs values
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:34
Data operator to do calculations at integration points.
Interface for managing meshsets containing materials and boundary conditions.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:665
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:315
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:387
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:639
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:269
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:713
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:490
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:482
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:287
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:459
std::vector< std::string > & getOtherFiniteElements()
Get the Other Finite Elements.
Definition: Simple.hpp:367
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:308
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:34
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Class implemented by user to detect face orientation.
Shape preserving constrains, i.e. nodes sliding on body surface.
Volume Length Quality.