30int main(
int argc,
char *argv[]) {
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,
48 CHKERR PetscOptionsScalar(
"-gamma_factor",
"",
51 ierr = PetscOptionsEnd();
56 moab::Interface &moab = moab_core;
107 ->synchroniseFieldEntities(
"LAMBDA_EDGE");
112 "EDGE_SLIDING",
"MESH_NODE_POSITIONS");
114 "EDGE_SLIDING",
"MESH_NODE_POSITIONS");
116 "EDGE_SLIDING",
"MESH_NODE_POSITIONS");
146 struct ElementsAndOperators {
150 double *minQualityPtr;
155 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
156 ierr = VecGetArray(minQualityVec, &minQualityPtr);
157 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
160 virtual ~ElementsAndOperators() {
161 ierr = VecRestoreArray(minQualityVec, &minQualityPtr);
162 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
163 ierr = VecDestroy(&minQualityVec);
164 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
167 double getMinQuality()
const {
return *minQualityPtr; }
171 SURFACE_CONSTRAINS_TAG,
175 boost::shared_ptr<Smoother> smootherFe;
176 boost::shared_ptr<Smoother::MyVolumeFE>
178 boost::shared_ptr<Smoother::MyVolumeFE>
180 boost::shared_ptr<VolumeLengthQuality<double> > volumeLengthDouble;
181 boost::shared_ptr<VolumeLengthQuality<adouble> > volumeLengthAdouble;
183 boost::shared_ptr<SurfaceSlidingConstrains> surfaceConstrain;
184 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
186 boost::shared_ptr<EdgeSlidingConstrains> edgeConstrain;
188 boost::shared_ptr<DirichletFixFieldAtEntitiesBc> fixMaterialEnts;
190 boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore> minQualityFe;
193 double *minQualityPtr;
194 MinQualityOp(
double *min_quality_ptr)
196 "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
197 minQualityPtr(min_quality_ptr) {}
201 if (type != MBVERTEX)
204 *minQualityPtr = fmin(*minQualityPtr, q);
211 smootherFe = boost::shared_ptr<Smoother>(
new Smoother(mField));
212 volumeLengthAdouble = boost::shared_ptr<VolumeLengthQuality<adouble> >(
214 volumeLengthDouble = boost::shared_ptr<VolumeLengthQuality<double> >(
219 smootherFe->setOfBlocks[0].tEts.merge(tets);
221 smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
222 smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
225 smootherFe->commonData.spatialPositions =
"MESH_NODE_POSITIONS";
226 smootherFe->commonData.meshPositions =
"NONE";
228 smootherFe->feRhs.meshPositionsFieldName =
"NONE";
229 smootherFe->feLhs.meshPositionsFieldName =
"NONE";
230 smootherFe->feRhs.addToRule = 0;
231 smootherFe->feLhs.addToRule = 0;
233 feSmootherRhs = smootherFe->feRhsPtr;
234 feSmootherLhs = smootherFe->feLhsPtr;
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));
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"));
262 boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore>(
264 minQualityFe->getOpPtrVector().push_back(
265 new MinQualityOp(minQualityPtr));
273 fixMaterialEnts = boost::shared_ptr<DirichletFixFieldAtEntitiesBc>(
276 fixMaterialEnts->fieldNames.push_back(
"LAMBDA_SURFACE");
277 fixMaterialEnts->fieldNames.push_back(
"LAMBDA_EDGE");
284 skinOrientation = boost::shared_ptr<
287 surfaceConstrain = boost::shared_ptr<SurfaceSlidingConstrains>(
290 surfaceConstrain->setOperators(SURFACE_CONSTRAINS_TAG,
"LAMBDA_SURFACE",
291 "MESH_NODE_POSITIONS");
304 CHKERR skin.find_skin(0, tets,
false, skin_faces);
306 edgeConstrain = boost::shared_ptr<EdgeSlidingConstrains>(
308 CHKERR edgeConstrain->setOperators(EDGE_CONSTRAINS_TAG, edges,
309 skin_faces,
"LAMBDA_EDGE",
310 "MESH_NODE_POSITIONS");
321 boost::shared_ptr<ForcesAndSourcesCore> null;
328 surfaceConstrain->feRhsPtr, null, null);
330 edgeConstrain->feRhsPtr, null, null);
339 surfaceConstrain->feLhsPtr, null, null);
341 edgeConstrain->feLhsPtr, null, null);
356 CHKERR VecMin(minQualityVec, PETSC_NULL, minQualityPtr);
361 ElementsAndOperators elements_and_operators(m_field);
362 CHKERR elements_and_operators.createSmoothingFE();
363 CHKERR elements_and_operators.createConstrians();
367 CHKERR elements_and_operators.addFEtoDM(dm);
376 CHKERR DMCreateGlobalVector(dm, &
F);
381 CHKERR zeroLambdaFields(dm);
386 CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
387 CHKERR SNESSetFromOptions(solver);
388 CHKERR SNESSetDM(solver, dm);
400 CHKERR SNESDestroy(&solver);
412 "MESH_NODE_POSITIONS", it)) {
413 if (it->get()->getEntType() != MBVERTEX)
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());
429 "MESH_NODE_POSITIONS", it)) {
430 if (it->get()->getEntType() != MBVERTEX)
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];
447 0, MBVERTEX,
"LAMBDA_SURFACE");
454 CHKERR solve.setFieldFromCoords(dm);
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);
460 double gamma = min_quality > 0 ?
gamma_factor * min_quality
462 elements_and_operators.volumeLengthDouble->gAmma = gamma;
463 elements_and_operators.volumeLengthAdouble->gAmma = gamma;
465 double min_quality_p,
eps;
468 min_quality_p = min_quality;
472 CHKERR solve.setCoordsFromField(dm);
473 CHKERR elements_and_operators.calcuteMinQuality(dm);
474 min_quality = elements_and_operators.getMinQuality();
476 eps = (min_quality - min_quality_p) / min_quality;
477 PetscPrintf(PETSC_COMM_WORLD,
"Min quality = %4.3f eps = %4.3f\n",
480 double gamma = min_quality > 0 ?
gamma_factor * min_quality
482 elements_and_operators.volumeLengthDouble->gAmma = gamma;
483 elements_and_operators.volumeLengthAdouble->gAmma = gamma;