42 {
44
45 try {
46
47 moab::Core mb_instance;
48 moab::Interface& moab = mb_instance;
49
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;
56
57
58 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Bone remodeling",
"none"); CHKERRQ(
ierr);
59 ierr = PetscOptionsString(
60 "-my_file",
61 "mesh file name","",
64
65 ierr = PetscOptionsString(
66 "-meta_image",
67 "meta image name","",
68 "file.mhd",meta_file_name,255,&meta_flg_file
70
71
72
73
74
75
76
77 ierr = PetscOptionsBool(
78 "-my_is_atom_test",
79 "is used with testing, exit with error when diverged","",
80 PETSC_FALSE,&is_atom_test,PETSC_NULL
82
83 ierr = PetscOptionsInt(
84 "-my_order",
85 "default approximation order","",
88
89 ierr = PetscOptionsEnd(); CHKERRQ(
ierr);
90
91
92
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] << " ";
106 }
107 cout << endl;
108 }
109
110
111 if(flg_file != PETSC_TRUE) {
112 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
113 }
114 const char *option;
115 option = "";
117 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
118 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
119
122
124
125
127 bit_level0.set(0);
129
132
133
135 ierr = meshsets_manager_ptr->setMeshsetFromFile(); CHKERRQ(
ierr);
136
137 PetscPrintf(PETSC_COMM_WORLD,"Read meshset. Added meshsets for bc.cfg\n");
139 PetscPrintf(
140 PETSC_COMM_WORLD,
141 "%s",static_cast<std::ostringstream&>(std::ostringstream().seekp(0) << *it << endl).str().c_str()
142 );
143 cerr << *it << endl;
144 }
145
148
150 string name = it->getName();
151 if (name.compare(0,14,"NO_REMODELLING") == 0) {
152
156 meshset,MBTET,tets_block,true
158 tets_all = subtract(tets_all,tets_block);
159 }
160 }
161
162
163
164
165
166
167
168
169
170
171
172
175
176
177 {
187 }
188
189
196
197
198 {
199
200 Skinner skin(&moab);
205 rval = moab.get_adjacencies(
206 faces,1,false,edges,moab::Interface::UNION
209 }
211
213
214 {
216 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(
ierr);
217 }
218
220 {
223 ierr = surface_distance.takeASkin(tets); CHKERRQ(
ierr);
224 ierr = surface_distance.buildTree(); CHKERRQ(
ierr);
225 ierr = surface_distance.setDistanceFromSurface(
"DISTANCE_FROM_SURFACE"); CHKERRQ(
ierr);
226 }
227
228
235
236
238
239
241
243
244 DM dm_rho;
245 DMType dm_name = "DM_RHO";
246
248 ierr = DMCreate(PETSC_COMM_WORLD,&dm_rho);CHKERRQ(
ierr);
249 ierr = DMSetType(dm_rho,dm_name);CHKERRQ(
ierr);
250
252
253 ierr = DMSetFromOptions(dm_rho); CHKERRQ(
ierr);
254
256
257
258
259 {
260 ierr = DMSetUp(dm_rho); CHKERRQ(
ierr);
261
262
263
264
265
266
267
268
269 }
270
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);
281
283 auto &list_of_operators = fe_volume.getOpPtrVector();
286 ierr = data_from_meta_io.fromOptions(); CHKERRQ(
ierr);
287 list_of_operators.push_back(
new OpMassCalculation(
"RHO",mass_vec,data_from_meta_io));
289
290 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(
ierr);
291 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(
ierr);
292 double volume;
293 ierr = VecSum(volume_vec,&volume); CHKERRQ(
ierr);
294 ierr = PetscPrintf(PETSC_COMM_WORLD,
"Volume %0.32g\n",volume); CHKERRQ(
ierr);
295
296 ierr = VecAssemblyBegin(mass_vec); CHKERRQ(
ierr);
297 ierr = VecAssemblyEnd(mass_vec); CHKERRQ(
ierr);
298 double mass;
299 ierr = VecSum(mass_vec,&mass); CHKERRQ(
ierr);
300 ierr = PetscPrintf(PETSC_COMM_WORLD,
"Mass %0.32g\n",mass); CHKERRQ(
ierr);
301
302 ierr = VecDestroy(&volume_vec); CHKERRQ(
ierr);
303 ierr = VecDestroy(&mass_vec); CHKERRQ(
ierr);
304
305
308 ierr = DMCreateMatrix(dm_rho,&A); CHKERRQ(
ierr);
311
312 ierr = MatZeroEntries(A); CHKERRQ(
ierr);
314
316 common_data.globA =
A;
317 common_data.globF =
F;
318
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;
323
325
326 fe_density_approximation.getOpPtrVector().push_back
328 fe_density_approximation.getOpPtrVector().push_back
330 fe_density_approximation.getOpPtrVector().push_back
332 fe_density_approximation.getOpPtrVector().push_back
334 fe_density_approximation.getOpPtrVector().push_back
337
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);
344
345
346
347
348
349 {
350 KSP ksp;
351 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(
ierr);
352 ierr = KSPSetOperators(ksp,A,A); CHKERRQ(
ierr);
353 ierr = KSPSetFromOptions(ksp); CHKERRQ(
ierr);
355 }
356
357 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
358 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
360
362
363 fe_post_process_calulations.getOpPtrVector().push_back(
365 );
367
368 ierr = VecAssemblyBegin(mass_vec_approx); CHKERRQ(
ierr);
369 ierr = VecAssemblyEnd(mass_vec_approx); CHKERRQ(
ierr);
370 double mass_approx;
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);
373 if(is_atom_test) {
374 if( abs(mass_approx - volume) > 1.0e-2) {
376 }
377 }
378 ierr = VecDestroy(&mass_vec_approx); CHKERRQ(
ierr);
379
381 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(
ierr);
382 ierr = post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS"); CHKERRQ(
ierr);
383 ierr = post_proc.addFieldValuesPostProc(
"RHO"); CHKERRQ(
ierr);
384 ierr = post_proc.addFieldValuesPostProc(
"DISTANCE_FROM_SURFACE"); CHKERRQ(
ierr);
385 ierr = post_proc.addFieldValuesGradientPostProc(
"RHO"); CHKERRQ(
ierr);
386 ierr = post_proc.addFieldValuesGradientPostProc(
"DISTANCE_FROM_SURFACE"); CHKERRQ(
ierr);
388 ierr = post_proc.writeFile(
"out.h5m"); CHKERRQ(
ierr);
389
392 ierr = MatDestroy(&A); CHKERRQ(
ierr);
393 ierr = DMDestroy(&dm_rho);
394
395
398 CHKERR moab.write_file(
"analysis_mesh.h5m");
399
400
402
404
405 return (0);
406}
#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.
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.
const FTensor::Tensor2< T, Dim, Dim > Vec
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
Assemble local vector containing density data.
Calculate mass after approximation.
Calculate mass before approximation.
Create KDTree on surface of the mesh and calculate distance.
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...
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.