#include <Smoother.hpp>
#include <SurfaceSlidingConstrains.hpp>
#include <VolumeLengthQuality.hpp>
static char help[] =
"mesh smoothing\n\n";
int main(
int argc,
char *argv[]) {
try {
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Mesh cut options",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsInt(
"-edges_block_set",
"edges side set",
"",
CHKERR PetscOptionsInt(
"-vertex_block_set",
"vertex side set",
"",
CHKERR PetscOptionsBool(
"-output_vtk",
"if true outout vtk file",
"",
CHKERR PetscOptionsScalar(
"-quality_reduction_tol",
"",
"Tolerance of quality reduction",
tol, &
tol,
PETSC_NULL);
CHKERR PetscOptionsScalar(
"-gamma_factor",
"",
PETSC_NULL);
ierr = PetscOptionsEnd();
{
1);
1);
{
true);
"LAMBDA_EDGE");
->synchroniseFieldEntities("LAMBDA_EDGE");
"EDGE_SLIDING", "MESH_NODE_POSITIONS");
"EDGE_SLIDING", "MESH_NODE_POSITIONS");
"EDGE_SLIDING", "MESH_NODE_POSITIONS");
"LAMBDA_EDGE");
"LAMBDA_EDGE");
"LAMBDA_EDGE");
"EDGE_SLIDING");
}
}
}
}
struct ElementsAndOperators {
double *minQualityPtr;
&minQualityVec);
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
ierr = VecGetArray(minQualityVec, &minQualityPtr);
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
virtual ~ElementsAndOperators() {
ierr = VecRestoreArray(minQualityVec, &minQualityPtr);
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
ierr = VecDestroy(&minQualityVec);
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
double getMinQuality() const { return *minQualityPtr; }
enum Tags {
SURFACE_CONSTRAINS_TAG,
EDGE_CONSTRAINS_TAG
};
boost::shared_ptr<Smoother> smootherFe;
boost::shared_ptr<Smoother::MyVolumeFE>
feSmootherRhs;
boost::shared_ptr<Smoother::MyVolumeFE>
feSmootherLhs;
boost::shared_ptr<VolumeLengthQuality<double> > volumeLengthDouble;
boost::shared_ptr<VolumeLengthQuality<adouble> > volumeLengthAdouble;
boost::shared_ptr<SurfaceSlidingConstrains> surfaceConstrain;
boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
skinOrientation;
boost::shared_ptr<EdgeSlidingConstrains> edgeConstrain;
boost::shared_ptr<DirichletFixFieldAtEntitiesBc> fixMaterialEnts;
boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore> minQualityFe;
double *minQualityPtr;
MinQualityOp(double *min_quality_ptr)
minQualityPtr(min_quality_ptr) {}
*minQualityPtr = fmin(*minQualityPtr, q);
}
};
smootherFe = boost::shared_ptr<Smoother>(
new Smoother(mField));
volumeLengthAdouble = boost::shared_ptr<VolumeLengthQuality<adouble> >(
volumeLengthDouble = boost::shared_ptr<VolumeLengthQuality<double> >(
smootherFe->setOfBlocks[0].tEts.merge(tets);
smootherFe->setOfBlocks[0].materialDoublePtr = volumeLengthDouble;
smootherFe->setOfBlocks[0].materialAdoublePtr = volumeLengthAdouble;
smootherFe->commonData.spatialPositions = "MESH_NODE_POSITIONS";
smootherFe->commonData.meshPositions = "NONE";
smootherFe->feRhs.meshPositionsFieldName = "NONE";
smootherFe->feLhs.meshPositionsFieldName = "NONE";
smootherFe->feRhs.addToRule = 0;
smootherFe->feLhs.addToRule = 0;
feSmootherRhs = smootherFe->feRhsPtr;
feSmootherLhs = smootherFe->feLhsPtr;
feSmootherRhs->getOpPtrVector().push_back(
"MESH_NODE_POSITIONS", smootherFe->commonData));
feSmootherRhs->getOpPtrVector().push_back(
"MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
"MESH_NODE_POSITIONS", smootherFe->setOfBlocks[0],
smootherFe->commonData, smootherFe->smootherData));
feSmootherLhs->getOpPtrVector().push_back(
"MESH_NODE_POSITIONS", smootherFe->commonData));
feSmootherLhs->getOpPtrVector().push_back(
"MESH_NODE_POSITIONS", smootherFe->setOfBlocks.at(0),
"MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
smootherFe->setOfBlocks.at(0), smootherFe->commonData,
smootherFe->smootherData, "LAMBDA_CRACKFRONT_AREA_TANGENT"));
minQualityFe =
boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore>(
minQualityFe->getOpPtrVector().push_back(
new MinQualityOp(minQualityPtr));
}
fixMaterialEnts = boost::shared_ptr<DirichletFixFieldAtEntitiesBc>(
fixed_vertex));
fixMaterialEnts->fieldNames.push_back("LAMBDA_SURFACE");
fixMaterialEnts->fieldNames.push_back("LAMBDA_EDGE");
}
skinOrientation = boost::shared_ptr<
surfaceConstrain = boost::shared_ptr<SurfaceSlidingConstrains>(
skinOrientation,
surfaceConstrain->setOperators(SURFACE_CONSTRAINS_TAG, "LAMBDA_SURFACE",
"MESH_NODE_POSITIONS");
true);
CHKERR skin.find_skin(0, tets,
false, skin_faces);
edgeConstrain = boost::shared_ptr<EdgeSlidingConstrains>(
CHKERR edgeConstrain->setOperators(EDGE_CONSTRAINS_TAG, edges,
skin_faces, "LAMBDA_EDGE",
"MESH_NODE_POSITIONS");
}
}
boost::shared_ptr<ForcesAndSourcesCore> null;
null);
null);
surfaceConstrain->feRhsPtr, null, null);
edgeConstrain->feRhsPtr, null, null);
fixMaterialEnts);
null);
null);
surfaceConstrain->feLhsPtr, null, null);
edgeConstrain->feLhsPtr, null, null);
fixMaterialEnts);
}
*minQualityPtr = 1;
CHKERR VecMin(minQualityVec, PETSC_NULL, minQualityPtr);
}
};
ElementsAndOperators elements_and_operators(m_field);
CHKERR elements_and_operators.createSmoothingFE();
CHKERR elements_and_operators.createConstrians();
DM dm;
CHKERR elements_and_operators.addFEtoDM(dm);
struct Solve {
CHKERR DMCreateGlobalVector(dm, &
F);
SNES solver;
CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
CHKERR SNESSetFromOptions(solver);
}
"MESH_NODE_POSITIONS", it)) {
if (it->get()->getEntType() != MBVERTEX)
continue;
coords[
dd] = it->get()->getEntFieldData()[
dd];
}
}
"MESH_NODE_POSITIONS", it)) {
if (it->get()->getEntType() != MBVERTEX)
continue;
it->get()->getEntFieldData()[
dd] = coords[
dd];
}
}
private:
0, MBVERTEX, "LAMBDA_SURFACE");
}
};
CHKERR elements_and_operators.calcuteMinQuality(dm);
double min_quality = elements_and_operators.getMinQuality();
PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f\n", min_quality);
elements_and_operators.volumeLengthDouble->gAmma = gamma;
elements_and_operators.volumeLengthAdouble->gAmma = gamma;
double min_quality_p,
eps;
do {
min_quality_p = min_quality;
CHKERR elements_and_operators.calcuteMinQuality(dm);
min_quality = elements_and_operators.getMinQuality();
eps = (min_quality - min_quality_p) / min_quality;
PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f eps = %4.3f\n",
elements_and_operators.volumeLengthDouble->gAmma = gamma;
elements_and_operators.volumeLengthAdouble->gAmma = gamma;
"");
}
}