![Logo](MoFEMLogo.png) |
| v0.14.0
|
Go to the documentation of this file.
10 using namespace MoFEM;
18 using namespace MoFEM;
20 static char help[] =
"mesh smoothing\n\n";
30 int 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();
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)
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),
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),
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");
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;
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;
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Data on single entity (This is passed as argument to DataOperator::doWork)
VectorBoundedArray< double, 3 > VectorDouble3
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
MoFEMErrorCode buildProblem()
Build problem.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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
const VectorDouble & getFieldData() const
get dofs values
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Simple interface for fast problem set-up.
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.
Deprecated interface functions.
DeprecatedCoreInterface Interface
MoFEMErrorCode getOptions()
get options
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
MoFEMErrorCode getDM(DM *dm)
Get DM.
#define CHKERR
Inline error check.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual moab::Interface & get_moab()=0
implementation of Data Operators for Forces and Sources
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_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.
int main(int argc, char *argv[])
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
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.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
virtual int get_comm_size() const =0
const std::string getDomainFEName() const
Get the Domain FE Name.
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Shape preserving constrains, i.e. nodes sliding on body surface.
friend class UserDataOperator
Volume finite element base.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode buildFiniteElements()
Build finite elements.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
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)
#define CATCH_ERRORS
Catch errors.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
const FTensor::Tensor2< T, Dim, Dim > Vec
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode defineFiniteElements()
Define finite elements.
MoFEMErrorCode buildFields()
Build fields.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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.
Implementation of Volume-Length-Quality measure with barrier.
const double D
diffusivity
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Implementing surface sliding constrains.
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.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
std::vector< std::string > & getOtherFiniteElements()
Get the Other Finite Elements.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Class implemented by user to detect face orientation.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ OPROW
operator doWork function is executed on FE rows
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function