42 {
43
45
46 try {
47
48 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Mesh cut options",
"none");
49 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
51 CHKERR PetscOptionsInt(
"-edges_block_set",
"edges side set",
"",
53 CHKERR PetscOptionsInt(
"-vertex_block_set",
"vertex side set",
"",
55 CHKERR PetscOptionsBool(
"-output_vtk",
"if true outout vtk file",
"",
57 CHKERR PetscOptionsScalar(
"-quality_reduction_tol",
"",
58 "Tolerance of quality reduction",
tol, &
tol,
59 PETSC_NULL);
60 CHKERR PetscOptionsScalar(
"-gamma_factor",
"",
62 PETSC_NULL);
63 ierr = PetscOptionsEnd();
65
66
69
72
74
75
76
78
80
81
82 {
83
85
87
88
89
92
95
96
98 1);
100 1);
101
104
105
106 {
109
112 Range edges;
115 true);
117 "LAMBDA_EDGE");
119 ->synchroniseFieldEntities("LAMBDA_EDGE");
121
124 "EDGE_SLIDING", "MESH_NODE_POSITIONS");
126 "EDGE_SLIDING", "MESH_NODE_POSITIONS");
128 "EDGE_SLIDING", "MESH_NODE_POSITIONS");
130 "LAMBDA_EDGE");
132 "LAMBDA_EDGE");
134 "LAMBDA_EDGE");
136 "EDGE_SLIDING");
138 }
139 }
140
144
147 Range edges;
150 Range verts;
153 }
156 }
157
158 struct ElementsAndOperators {
159
162 double *minQualityPtr;
163
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 {
183 SURFACE_CONSTRAINS_TAG,
184 EDGE_CONSTRAINS_TAG
185 };
186
187 boost::shared_ptr<Smoother> smootherFe;
188 boost::shared_ptr<Smoother::MyVolumeFE>
189 feSmootherRhs;
190 boost::shared_ptr<Smoother::MyVolumeFE>
191 feSmootherLhs;
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
205 double *minQualityPtr;
206 MinQualityOp(double *min_quality_ptr)
209 minQualityPtr(min_quality_ptr) {}
213 if (
type != MBVERTEX)
215 double q = Tools::volumeLengthQuality(&*data.
getFieldData().begin());
216 *minQualityPtr = fmin(*minQualityPtr, q);
218 }
219 };
220
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;
231 smootherFe->setOfBlocks[0].tEts.merge(tets);
232
233 smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
234 smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
235
236
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
249 feSmootherRhs->getOpPtrVector().push_back(
251 "MESH_NODE_POSITIONS", smootherFe->commonData));
252 feSmootherRhs->getOpPtrVector().push_back(
254 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
257 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks[0],
258 smootherFe->commonData, smootherFe->smootherData));
259
260
261 feSmootherLhs->getOpPtrVector().push_back(
263 "MESH_NODE_POSITIONS", smootherFe->commonData));
264 feSmootherLhs->getOpPtrVector().push_back(
266 "MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
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;
284 }
285 fixMaterialEnts = boost::shared_ptr<DirichletFixFieldAtEntitiesBc>(
287 fixed_vertex));
288 fixMaterialEnts->fieldNames.push_back("LAMBDA_SURFACE");
289 fixMaterialEnts->fieldNames.push_back("LAMBDA_EDGE");
290
292 }
293
296 skinOrientation = boost::shared_ptr<
299 surfaceConstrain = boost::shared_ptr<SurfaceSlidingConstrains>(
300 skinOrientation,
302 surfaceConstrain->setOperators(SURFACE_CONSTRAINS_TAG, "LAMBDA_SURFACE",
303 "MESH_NODE_POSITIONS");
304
307 Range edges;
310 true);
311
312 Range tets;
315 Range skin_faces;
316 CHKERR skin.find_skin(0, tets,
false, skin_faces);
317
318 edgeConstrain = boost::shared_ptr<EdgeSlidingConstrains>(
320 CHKERR edgeConstrain->setOperators(EDGE_CONSTRAINS_TAG, edges,
321 skin_faces, "LAMBDA_EDGE",
322 "MESH_NODE_POSITIONS");
323
324
325
326
327 }
329 }
330
333 boost::shared_ptr<ForcesAndSourcesCore> null;
334
336 null);
338 null);
340 surfaceConstrain->feRhsPtr, null, null);
342 edgeConstrain->feRhsPtr, null, null);
344 fixMaterialEnts);
345
347 null);
349 null);
351 surfaceConstrain->feLhsPtr, null, null);
353 edgeConstrain->feLhsPtr, null, null);
355 fixMaterialEnts);
356
357
358
359
360
362 }
363
366 *minQualityPtr = 1;
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;
379 CHKERR elements_and_operators.addFEtoDM(dm);
380
381 struct Solve {
382
385
386
388 CHKERR DMCreateGlobalVector(dm, &F);
389
390
392
393 CHKERR zeroLambdaFields(dm);
395
396
397 SNES solver;
398 CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
399 CHKERR SNESSetFromOptions(solver);
400 CHKERR SNESSetDM(solver, dm);
401
403
404
405
406 CHKERR SNESSolve(solver, F,
D);
407
408
409
411
412 CHKERR SNESDestroy(&solver);
415
417 }
418
424 "MESH_NODE_POSITIONS", it)) {
425 if (it->get()->getEntType() != MBVERTEX)
426 continue;
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
441 "MESH_NODE_POSITIONS", it)) {
442 if (it->get()->getEntType() != MBVERTEX)
443 continue;
444 EntityHandle ent = it->get()->getEnt();
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:
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
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
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",
491
492 double gamma = min_quality > 0 ?
gamma_factor * min_quality
494 elements_and_operators.volumeLengthDouble->gAmma = gamma;
495 elements_and_operators.volumeLengthAdouble->gAmma = gamma;
496
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
520 "");
521 }
523
525}
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 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
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
DeprecatedCoreInterface Interface
const double D
diffusivity
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.