#include <BasicFiniteElements.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();
Simple *simple_interface;
{
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile();
CHKERR simple_interface->addDomainField(
"MESH_NODE_POSITIONS",
H1,
CHKERR simple_interface->addBoundaryField(
"LAMBDA_SURFACE",
H1,
CHKERR simple_interface->setFieldOrder(
"MESH_NODE_POSITIONS",
1);
CHKERR simple_interface->setFieldOrder(
"LAMBDA_SURFACE",
1);
simple_interface->getDomainFEName() = "SMOOTHING";
simple_interface->getBoundaryFEName() = "SURFACE_SLIDING";
{
Range edges;
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");
simple_interface->getOtherFiniteElements().push_back("EDGE_SLIDING");
}
}
CHKERR simple_interface->defineFiniteElements();
CHKERR simple_interface->defineProblem();
CHKERR simple_interface->buildFields();
Range edges;
Range verts;
}
CHKERR simple_interface->buildFiniteElements();
CHKERR simple_interface->buildProblem();
}
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)
"MESH_NODE_POSITIONS", UserDataOperator::OPROW),
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> >(
Range tets;
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));
Range fixed_vertex;
}
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");
Range edges;
true);
Range tets;
Range skin_faces;
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 simple_interface->getDM(&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];
EntityHandle ent = it->get()->getEnt();
}
}
"MESH_NODE_POSITIONS", it)) {
if (it->get()->getEntType() != MBVERTEX)
continue;
EntityHandle ent = it->get()->getEnt();
it->get()->getEntFieldData()[
dd] = coords[
dd];
}
}
private:
0, MBVERTEX, "LAMBDA_SURFACE");
}
};
Solve solve;
CHKERR solve.setFieldFromCoords(dm);
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 solve.setCoordsFromField(dm);
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;
"");
}
}
Implementing surface sliding constrains.
Implementation of Volume-Length-Quality measure with barrier.
#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 DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
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 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
int main(int argc, char *argv[])
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.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
VectorBoundedArray< double, 3 > VectorDouble3
implementation of Data Operators for Forces and Sources
DeprecatedCoreInterface Interface
const double D
diffusivity
DataForcesAndSourcesCore::EntData EntData
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
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 operator to do calculations at integration points.
friend class UserDataOperator
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Volume finite element with switches.
Class implemented by user to detect face orientation.
Shape preserving constrains, i.e. nodes sliding on body surface.