25#include <moab/Skinner.hpp>
26#include <moab/AdaptiveKDTree.hpp>
39 #error "You need compile users modules with MetaIO -DWITH_METAIO=1"
42int main(
int argc,
char *argv[]) {
47 moab::Core mb_instance;
48 moab::Interface& moab = mb_instance;
51 PetscBool flg_file = PETSC_FALSE;
52 char meta_file_name[255];
53 PetscBool meta_flg_file = PETSC_FALSE;
55 PetscBool is_atom_test = PETSC_FALSE;
58 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling",
"none"); CHKERRQ(
ierr);
59 ierr = PetscOptionsString(
65 ierr = PetscOptionsString(
68 "file.mhd",meta_file_name,255,&meta_flg_file
77 ierr = PetscOptionsBool(
79 "is used with testing, exit with error when diverged",
"",
80 PETSC_FALSE,&is_atom_test,PETSC_NULL
83 ierr = PetscOptionsInt(
85 "default approximation order",
"",
89 ierr = PetscOptionsEnd(); CHKERRQ(
ierr);
93 MetaImage meta_file(meta_file_name);
94 cout <<
"Quantity " << meta_file.Quantity() << endl;
95 const int *dim_size = meta_file.DimSize();
96 cout <<
"Image dimension: " << dim_size[0] <<
" " << dim_size[1] <<
" " << dim_size[2] << endl;
97 const double *elemet_size = meta_file.ElementSize();
98 cout <<
"Image element size: " << elemet_size[0] <<
" " << elemet_size[1] <<
" " << elemet_size[2] << endl;
99 const double *offset = meta_file.Offset();
100 cout <<
"Image offset: " << offset[0] <<
" " << offset[1] <<
" " << offset[2] << endl;
101 const double * transform_matrix = meta_file.TransformMatrix();
102 cout <<
"Image Transform Matrix:\n";
103 for(
int rr = 0;rr<3;rr++) {
104 for(
int cc = 0;cc<3;cc++) {
105 cout << transform_matrix[rr*3+cc] <<
" ";
111 if(flg_file != PETSC_TRUE) {
112 SETERRQ(PETSC_COMM_SELF,1,
"*** ERROR -my_file (MESH FILE NEEDED)");
117 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
118 if(pcomm == NULL) pcomm =
new ParallelComm(&moab,PETSC_COMM_WORLD);
135 ierr = meshsets_manager_ptr->setMeshsetFromFile(); CHKERRQ(
ierr);
137 PetscPrintf(PETSC_COMM_WORLD,
"Read meshset. Added meshsets for bc.cfg\n");
141 "%s",
static_cast<std::ostringstream&
>(std::ostringstream().seekp(0) << *it << endl).str().c_str()
150 string name = it->getName();
151 if (name.compare(0,14,
"NO_REMODELLING") == 0) {
156 meshset,MBTET,tets_block,
true
158 tets_all = subtract(tets_all,tets_block);
205 rval = moab.get_adjacencies(
206 faces,1,
false,edges,moab::Interface::UNION
216 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(
ierr);
245 DMType dm_name =
"DM_RHO";
248 ierr = DMCreate(PETSC_COMM_WORLD,&dm_rho);CHKERRQ(
ierr);
249 ierr = DMSetType(dm_rho,dm_name);CHKERRQ(
ierr);
253 ierr = DMSetFromOptions(dm_rho); CHKERRQ(
ierr);
260 ierr = DMSetUp(dm_rho); CHKERRQ(
ierr);
272 int volume_vec_ghost[] = { 0 };
273 ierr = VecCreateGhost(
274 PETSC_COMM_WORLD,(!m_field.
get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
276 ierr = VecZeroEntries(volume_vec); CHKERRQ(
ierr);
278 ierr = VecDuplicate(volume_vec,&mass_vec); CHKERRQ(
ierr);
280 ierr = VecDuplicate(volume_vec,&mass_vec_approx); CHKERRQ(
ierr);
287 list_of_operators.push_back(
new OpMassCalculation(
"RHO",mass_vec,data_from_meta_io));
290 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(
ierr);
291 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(
ierr);
293 ierr = VecSum(volume_vec,&volume); CHKERRQ(
ierr);
294 ierr = PetscPrintf(PETSC_COMM_WORLD,
"Volume %0.32g\n",volume); CHKERRQ(
ierr);
296 ierr = VecAssemblyBegin(mass_vec); CHKERRQ(
ierr);
297 ierr = VecAssemblyEnd(mass_vec); CHKERRQ(
ierr);
299 ierr = VecSum(mass_vec,&mass); CHKERRQ(
ierr);
300 ierr = PetscPrintf(PETSC_COMM_WORLD,
"Mass %0.32g\n",mass); CHKERRQ(
ierr);
302 ierr = VecDestroy(&volume_vec); CHKERRQ(
ierr);
303 ierr = VecDestroy(&mass_vec); CHKERRQ(
ierr);
308 ierr = DMCreateMatrix(dm_rho,&
A); CHKERRQ(
ierr);
316 common_data.globA =
A;
317 common_data.globF =
F;
319 common_data.distanceValue = boost::shared_ptr<VectorDouble>(
new VectorDouble());
320 common_data.distanceGradient = boost::shared_ptr<MatrixDouble>(
new MatrixDouble());
321 boost::shared_ptr<VectorDouble> dist_value = common_data.distanceValue;
322 boost::shared_ptr<MatrixDouble> grad_value = common_data.distanceGradient;
338 ierr = MatAssemblyBegin(
A,MAT_FINAL_ASSEMBLY); CHKERRQ(
ierr);
339 ierr = MatAssemblyEnd(
A,MAT_FINAL_ASSEMBLY); CHKERRQ(
ierr);
340 ierr = VecGhostUpdateBegin(
F,ADD_VALUES,SCATTER_REVERSE);
341 ierr = VecGhostUpdateEnd(
F,ADD_VALUES,SCATTER_REVERSE);
342 ierr = VecAssemblyBegin(
F); CHKERRQ(
ierr);
351 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(
ierr);
352 ierr = KSPSetOperators(ksp,
A,
A); CHKERRQ(
ierr);
353 ierr = KSPSetFromOptions(ksp); CHKERRQ(
ierr);
357 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
358 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
368 ierr = VecAssemblyBegin(mass_vec_approx); CHKERRQ(
ierr);
369 ierr = VecAssemblyEnd(mass_vec_approx); CHKERRQ(
ierr);
371 ierr = VecSum(mass_vec_approx,&mass_approx); CHKERRQ(
ierr);
372 ierr = PetscPrintf(PETSC_COMM_WORLD,
"Mass after approximation %0.32g\n",mass_approx); CHKERRQ(
ierr);
374 if( abs(mass_approx - volume) > 1.0e-2) {
378 ierr = VecDestroy(&mass_vec_approx); CHKERRQ(
ierr);
393 ierr = DMDestroy(&dm_rho);
398 CHKERR moab.write_file(
"analysis_mesh.h5m");
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
Post-process fields on refined mesh.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRQ_MOAB(a)
check error code of MoAB function
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
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.
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 build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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_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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
Assemble local vector containing density data.
Calculate mass after approximation.
Calculate mass before approximation.
Create KDTree on surface of the mesh and calculate distance.
PetscErrorCode setDistanceFromSurface(const std::string field_name)
Set distance to vertices on mesh.
PetscErrorCode takeASkin(Range &volume)
Take a skin from a mesh.
PetscErrorCode buildTree()
Build tree.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode delete_finite_element(const std::string name, int verb=DEFAULT_VERBOSITY)=0
delete finite element from mofem database
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 int get_comm_rank() const =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.
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.