30 {
31
33
34 try {
35
36 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Mesh cut options",
"none");
37 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
39 CHKERR PetscOptionsInt(
"-edges_block_set",
"edges side set",
"",
41 CHKERR PetscOptionsInt(
"-vertex_block_set",
"vertex side set",
"",
43 CHKERR PetscOptionsBool(
"-output_vtk",
"if true outout vtk file",
"",
45 CHKERR PetscOptionsScalar(
"-quality_reduction_tol",
"",
46 "Tolerance of quality reduction",
tol, &
tol,
47 PETSC_NULL);
48 CHKERR PetscOptionsScalar(
"-gamma_factor",
"",
50 PETSC_NULL);
51 ierr = PetscOptionsEnd();
53
54
55 moab::Core moab_core;
56 moab::Interface &moab = moab_core;
57
60
62
63
64
66
68
69
70 {
71
73
75
76
77
80
83
84
86 1);
88 1);
89
92
93
94 {
97
103 true);
105 "LAMBDA_EDGE");
107 ->synchroniseFieldEntities("LAMBDA_EDGE");
109
112 "EDGE_SLIDING", "MESH_NODE_POSITIONS");
114 "EDGE_SLIDING", "MESH_NODE_POSITIONS");
116 "EDGE_SLIDING", "MESH_NODE_POSITIONS");
118 "LAMBDA_EDGE");
120 "LAMBDA_EDGE");
122 "LAMBDA_EDGE");
124 "EDGE_SLIDING");
126 }
127 }
128
132
141 }
144 }
145
146 struct ElementsAndOperators {
147
150 double *minQualityPtr;
151
154 &minQualityVec);
155 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
156 ierr = VecGetArray(minQualityVec, &minQualityPtr);
157 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
158 }
159
160 virtual ~ElementsAndOperators() {
161 ierr = VecRestoreArray(minQualityVec, &minQualityPtr);
162 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
163 ierr = VecDestroy(&minQualityVec);
164 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
165 }
166
167 double getMinQuality() const { return *minQualityPtr; }
168
169 enum Tags {
171 SURFACE_CONSTRAINS_TAG,
172 EDGE_CONSTRAINS_TAG
173 };
174
175 boost::shared_ptr<Smoother> smootherFe;
176 boost::shared_ptr<Smoother::MyVolumeFE>
177 feSmootherRhs;
178 boost::shared_ptr<Smoother::MyVolumeFE>
179 feSmootherLhs;
180 boost::shared_ptr<VolumeLengthQuality<double> > volumeLengthDouble;
181 boost::shared_ptr<VolumeLengthQuality<adouble> > volumeLengthAdouble;
182
183 boost::shared_ptr<SurfaceSlidingConstrains> surfaceConstrain;
184 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
185 skinOrientation;
186 boost::shared_ptr<EdgeSlidingConstrains> edgeConstrain;
187
188 boost::shared_ptr<DirichletFixFieldAtEntitiesBc> fixMaterialEnts;
189
190 boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore> minQualityFe;
191
193 double *minQualityPtr;
194 MinQualityOp(double *min_quality_ptr)
197 minQualityPtr(min_quality_ptr) {}
201 if (type != MBVERTEX)
204 *minQualityPtr = fmin(*minQualityPtr, q);
206 }
207 };
208
211 smootherFe = boost::shared_ptr<Smoother>(
new Smoother(mField));
212 volumeLengthAdouble = boost::shared_ptr<VolumeLengthQuality<adouble> >(
214 volumeLengthDouble = boost::shared_ptr<VolumeLengthQuality<double> >(
216
219 smootherFe->setOfBlocks[0].tEts.merge(tets);
220
221 smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
222 smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
223
224
225 smootherFe->commonData.spatialPositions = "MESH_NODE_POSITIONS";
226 smootherFe->commonData.meshPositions = "NONE";
227
228 smootherFe->feRhs.meshPositionsFieldName = "NONE";
229 smootherFe->feLhs.meshPositionsFieldName = "NONE";
230 smootherFe->feRhs.addToRule = 0;
231 smootherFe->feLhs.addToRule = 0;
232
233 feSmootherRhs = smootherFe->feRhsPtr;
234 feSmootherLhs = smootherFe->feLhsPtr;
235
236
237 feSmootherRhs->getOpPtrVector().push_back(
239 "MESH_NODE_POSITIONS", smootherFe->commonData));
240 feSmootherRhs->getOpPtrVector().push_back(
242 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
243 smootherFe->commonData, SMOOTHING_TAG, false));
245 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks[0],
246 smootherFe->commonData, smootherFe->smootherData));
247
248
249 feSmootherLhs->getOpPtrVector().push_back(
251 "MESH_NODE_POSITIONS", smootherFe->commonData));
252 feSmootherLhs->getOpPtrVector().push_back(
254 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
255 smootherFe->commonData, SMOOTHING_TAG, true));
257 "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
258 smootherFe->setOfBlocks.at(0), smootherFe->commonData,
259 smootherFe->smootherData, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
260
261 minQualityFe =
262 boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore>(
264 minQualityFe->getOpPtrVector().push_back(
265 new MinQualityOp(minQualityPtr));
266
272 }
273 fixMaterialEnts = boost::shared_ptr<DirichletFixFieldAtEntitiesBc>(
275 fixed_vertex));
276 fixMaterialEnts->fieldNames.push_back("LAMBDA_SURFACE");
277 fixMaterialEnts->fieldNames.push_back("LAMBDA_EDGE");
278
280 }
281
284 skinOrientation = boost::shared_ptr<
287 surfaceConstrain = boost::shared_ptr<SurfaceSlidingConstrains>(
288 skinOrientation,
290 surfaceConstrain->setOperators(SURFACE_CONSTRAINS_TAG, "LAMBDA_SURFACE",
291 "MESH_NODE_POSITIONS");
292
298 true);
299
304 CHKERR skin.find_skin(0, tets,
false, skin_faces);
305
306 edgeConstrain = boost::shared_ptr<EdgeSlidingConstrains>(
308 CHKERR edgeConstrain->setOperators(EDGE_CONSTRAINS_TAG, edges,
309 skin_faces, "LAMBDA_EDGE",
310 "MESH_NODE_POSITIONS");
311
312
313
314
315 }
317 }
318
321 boost::shared_ptr<ForcesAndSourcesCore> null;
322
324 null);
326 null);
328 surfaceConstrain->feRhsPtr, null, null);
330 edgeConstrain->feRhsPtr, null, null);
332 fixMaterialEnts);
333
335 null);
337 null);
339 surfaceConstrain->feLhsPtr, null, null);
341 edgeConstrain->feLhsPtr, null, null);
343 fixMaterialEnts);
344
345
346
347
348
350 }
351
354 *minQualityPtr = 1;
356 CHKERR VecMin(minQualityVec, PETSC_NULL, minQualityPtr);
358 }
359 };
360
361 ElementsAndOperators elements_and_operators(m_field);
362 CHKERR elements_and_operators.createSmoothingFE();
363 CHKERR elements_and_operators.createConstrians();
364
365 DM dm;
367 CHKERR elements_and_operators.addFEtoDM(dm);
368
369 struct Solve {
370
373
374
376 CHKERR DMCreateGlobalVector(dm, &
F);
377
378
380
381 CHKERR zeroLambdaFields(dm);
383
384
385 SNES solver;
386 CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
387 CHKERR SNESSetFromOptions(solver);
388 CHKERR SNESSetDM(solver, dm);
389
391
392
393
395
396
397
399
400 CHKERR SNESDestroy(&solver);
403
405 }
406
412 "MESH_NODE_POSITIONS", it)) {
413 if (it->get()->getEntType() != MBVERTEX)
414 continue;
416 for(
int dd = 0;
dd!=3;++
dd)
417 coords[dd] = it->get()->getEntFieldData()[
dd];
419 CHKERR m_field_ptr->
get_moab().set_coords(&ent, 1, &*coords.begin());
420 }
422 }
423
429 "MESH_NODE_POSITIONS", it)) {
430 if (it->get()->getEntType() != MBVERTEX)
431 continue;
434 CHKERR m_field_ptr->
get_moab().get_coords(&ent, 1, &*coords.begin());
435 for(
int dd = 0;
dd!=3;++
dd)
436 it->get()->getEntFieldData()[
dd] = coords[
dd];
437 }
439 }
440
441 private:
447 0, MBVERTEX, "LAMBDA_SURFACE");
449 }
450
451 };
452
453 Solve solve;
454 CHKERR solve.setFieldFromCoords(dm);
455
456 CHKERR elements_and_operators.calcuteMinQuality(dm);
457 double min_quality = elements_and_operators.getMinQuality();
458 PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f\n", min_quality);
459
460 double gamma = min_quality > 0 ?
gamma_factor * min_quality
462 elements_and_operators.volumeLengthDouble->gAmma = gamma;
463 elements_and_operators.volumeLengthAdouble->gAmma = gamma;
464
465 double min_quality_p,
eps;
466 do {
467
468 min_quality_p = min_quality;
469
471
472 CHKERR solve.setCoordsFromField(dm);
473 CHKERR elements_and_operators.calcuteMinQuality(dm);
474 min_quality = elements_and_operators.getMinQuality();
475
476 eps = (min_quality - min_quality_p) / min_quality;
477 PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f eps = %4.3f\n",
479
480 double gamma = min_quality > 0 ?
gamma_factor * min_quality
482 elements_and_operators.volumeLengthDouble->gAmma = gamma;
483 elements_and_operators.volumeLengthAdouble->gAmma = gamma;
484
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
508 "");
509 }
511
513}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
#define _IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
loop over all dofs from a moFEM field and particular field
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 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 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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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
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)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorDouble & getFieldData() const
get dofs values
structure to get information form mofem into EntitiesFieldData
Interface for managing meshsets containing materials and boundary conditions.
Simple interface for fast problem set-up.
MoFEMErrorCode buildProblem()
Build problem.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
MoFEMErrorCode buildFiniteElements()
Build finite elements.
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.
std::vector< std::string > & getOtherFiniteElements()
Get the Other Finite Elements.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode buildFields()
Build fields.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
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.
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
const std::string getDomainFEName() const
Get the Domain FE Name.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
Class implemented by user to detect face orientation.
Shape preserving constrains, i.e. nodes sliding on body surface.