v0.14.0
Loading...
Searching...
No Matches
DensityMaps.hpp
Go to the documentation of this file.
1/** \file DensityMaps.hpp
2* \ingroup bone_remodelling
3
4* \brief Operator can be used with any volume element to calculate sum of
5* volumes of all volumes in the set
6
7*/
8
9/*
10* This file is part of MoFEM.
11* MoFEM is free software: you can redistribute it and/or modify it under
12* the terms of the GNU Lesser General Public License as published by the
13* Free Software Foundation, either version 3 of the License, or (at your
14* option) any later version.
15*
16* MoFEM is distributed in the hope that it will be useful, but WITHOUT
17* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19* License for more details.
20*
21* You should have received a copy of the GNU Lesser General Public
22* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23
24#ifndef __DENSITY_MAPS_HPP__
25#define __DENSITY_MAPS_HPP__
26
27namespace BoneRemodeling {
28 /**
29 * \brief Create KDTree on surface of the mesh and calculate distance.
30 * \ingroup bone_remodelling
31 */
33
36 AdaptiveKDTree kdTree;
37
38
40 mField(m_field),
41 kdTree(&m_field.get_moab()) {
42 }
43
44
45
46
47 Range sKin; //< Skin mesh, i.e. mesh on surface
48
49 /**
50 * \brief Take a skin from a mesh.
51 *
52 * @param volume Volume elements in the mesh
53 * @return Error code
54 */
55 PetscErrorCode takeASkin(Range &volume) {
56 PetscFunctionBegin;
57 Skinner skin(&mField.get_moab());
58 rval = skin.find_skin(0,volume,false,sKin); CHKERRQ_MOAB(rval);
59 // Only for debugging
60 // EntityHandle meshset;
61 // rval = mField.get_moab().create_meshset(MESHSET_SET,meshset);
62 // rval = mField.get_moab().add_entities(meshset,sKin); CHKERRQ(rval);
63 // rval = mField.get_moab().write_file("skin.vtk","VTK","",&meshset,1); CHKERRQ(rval);
64 PetscFunctionReturn(0);
65 }
66
67 /**
68 * \brief Build tree
69 *
70 * @return Error code
71 */
72 PetscErrorCode buildTree() {
73 PetscFunctionBegin;
74 rval = kdTree.build_tree(sKin,&kdTreeRootMeshset); CHKERRQ_MOAB(rval);
75 PetscFunctionReturn(0);
76 }
77
78 /**
79 * \brief Find point on surface which is closet to given point coordinates
80 * @param x 1st coordinate of point
81 * @param y 2nd coordinate of point
82 * @param z 3rd coordinate of point
83 * @param dx coordinate of distance vector
84 * @param dy coordinate of distance vector
85 * @param z coordinate of distance vector
86 * @return Error code
87 */
89 const double x,
90 const double y,
91 const double z,
92 double &dx,
93 double &dy,
94 double &dz
95 ) {
96 PetscFunctionBegin;
97 const double coords[] = {x,y,z};
98 double closest_point[3];
99 EntityHandle triangle;
100 rval = kdTree.closest_triangle(
102 coords,
103 closest_point,
104 triangle
105 ); CHKERRQ_MOAB(rval);
106 cblas_daxpy(3,-1,coords,1,closest_point,1);
107 dx = closest_point[0];
108 dy = closest_point[1];
109 dz = closest_point[2];
110 PetscFunctionReturn(0);
111 }
112
113 /**
114 * \brief Set distance to vertices on mesh
115 *
116 * @param field_name Name of field for distances
117 * @return Error Code
118 */
119 PetscErrorCode setDistanceFromSurface(const std::string field_name) {
120 PetscFunctionBegin;
121 auto field_ents = mField.get_field_ents();
122 auto field_bit_number = mField.get_field_bit_number(field_name);
123 auto it = field_ents->lower_bound(
124 FieldEntity::getLoBitNumberUId(field_bit_number));
125 auto hi_it = field_ents->upper_bound(
126 FieldEntity::getHiBitNumberUId(field_bit_number));
127 // Get values at nodes first
128 for(;it!=hi_it;it++) {
129 EntityType type = it->get()->getEntType();
130 if(type!=MBVERTEX) {
131 continue;
132 }
133 EntityHandle ent = it->get()->getEnt();
134 double coords[3];
135 rval = mField.get_moab().get_coords(&ent,1,coords); CHKERRQ_MOAB(rval);
136 double delta[3];
138 coords[0],coords[1],coords[2],delta[0],delta[1],delta[2]
139 ); CHKERRQ(ierr);
140 // Get vector of DOFs on vertex
141 VectorAdaptor dofs = it->get()->getEntFieldData();
142 // Set values
143 dofs[0] = cblas_dnrm2(3,delta,1);
144 }
145 // Get values at edges and project those values on approximation base
146 it = field_ents->lower_bound(
147 FieldEntity::getLoBitNumberUId(field_bit_number));
148 for(;it!=hi_it;it++) {
149 EntityType type = it->get()->getEntType();
150 if(type!=MBEDGE) {
151 continue;
152 }
153 const EntityHandle *conn;
154 int num_ndodes;
155 EntityHandle ent = it->get()->getEnt();
156 rval = mField.get_moab().get_connectivity(ent,conn,num_ndodes,true); CHKERRQ_MOAB(rval);
157 double coords[3*num_ndodes];
158 rval = mField.get_moab().get_coords(conn,num_ndodes,coords); CHKERRQ_MOAB(rval);
159 cblas_daxpy(3,1,&coords[3],1,coords,1);
160 cblas_dscal(3,0.5,coords,1);
161 double delta[3];
163 coords[0],coords[1],coords[2],delta[0],delta[1],delta[2]
164 ); CHKERRQ(ierr);
165
166 auto conn_it0 = field_ents->find(
167 FieldEntity::getLoLocalEntityBitNumber(field_bit_number, conn[0]));
168 auto conn_it1 = field_ents->find(
169 FieldEntity::getLoLocalEntityBitNumber(field_bit_number, conn[1]));
170
171 if(conn_it0==field_ents->end()) {
172 SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"entity not found");
173 }
174 if (conn_it1 == field_ents->end()) {
175 SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"entity not found");
176 }
177
178 VectorAdaptor vec_conn0 = conn_it0->get()->getEntFieldData();
179 VectorAdaptor vec_conn1 = conn_it1->get()->getEntFieldData();
180 // cerr << vec_conn0 << " " << vec_conn1 << endl;
181 double ave_delta = 0.5*(vec_conn0[0]+vec_conn1[0]);
182 double edge_shape_function_val = 0.25;
183 if(it->get()->getApproxBase()==AINSWORTH_LEGENDRE_BASE) {
184 } else if(it->get()->getApproxBase()==AINSWORTH_LOBATTO_BASE) {
185 edge_shape_function_val *= LOBATTO_PHI0(0);
186 } else {
187 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"Base not implemented");
188 }
189 // Get vector of DOFs on vertex
190 VectorAdaptor dofs = it->get()->getEntFieldData();
191 // cerr << mid_val_dx << " " << ave_dx << endl;
192 // Set values
193 dofs[0] = (cblas_dnrm2(3,delta,1)-ave_delta)/edge_shape_function_val;
194 }
195 PetscFunctionReturn(0);
196 }
197
198 };
199
200 #ifdef WITH_METAIO
201
202 /** \brief Load data from MetaImage file, translate grayscale values to densities.
203 * \ingroup bone_remodelling
204 */
206
207 MetaImage &metaFile;
210 }
211
212
213 double sIgma;
214 double cUbe_size;
216 double sCale;
217 double dIst_tol;
218 double tHreshold;
219 double paramElastic[2];
220 double paramDensity[2];
221 double lAmbda;
222
223 PetscErrorCode fromOptions() {
224 PetscFunctionBegin;
225 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Get parameters for DataFromMetaIO","none"); CHKERRQ(ierr);
226
227 PetscBool flg = PETSC_TRUE;
228
229 sIgma = 1;
230 ierr = PetscOptionsScalar(
231 "-density_map_sigma",
232 "get sigma for gauss weighted average","",
233 sIgma,&sIgma,PETSC_NULL
234 ); CHKERRQ(ierr);
235
236 cUbe_size = 0;
237 ierr = PetscOptionsScalar(
238 "-cube_size",
239 "get cube size for the averaging","",
240 0,&cUbe_size,PETSC_NULL
241 ); CHKERRQ(ierr);
242
243 dIst_tol = 0;
244 ierr = PetscOptionsScalar(
245 "-dist_tol",
246 "get distance tolerance from mesh surface","",
247 0,&dIst_tol,PETSC_NULL
248 ); CHKERRQ(ierr);
249
250 sCale = 1;
251 ierr = PetscOptionsScalar(
252 "-scale",
253 "scale the distance between voxels (eg. from mm to m)","",
254 sCale,&sCale,PETSC_NULL
255 ); CHKERRQ(ierr);
256
257 tHreshold = -1024;
258 ierr = PetscOptionsScalar(
259 "-threshold",
260 "minimum value of ct threshold","",
261 tHreshold,&tHreshold,PETSC_NULL
262 ); CHKERRQ(ierr);
263
264 fIlter = 0;
265 ierr = PetscOptionsInt(
266 "-filter",
267 "number of the filter","",
268 fIlter,&fIlter,PETSC_NULL
269 ); CHKERRQ(ierr);
270
271 lAmbda = 1;
272 ierr = PetscOptionsScalar(
273 "-lambda",
274 "get lambda coefficient for smoothing","",
275 lAmbda,&lAmbda,PETSC_NULL
276 ); CHKERRQ(ierr);
277
278 paramElastic[0] = 0;
279 paramElastic[1] = 1;
280 int nmax = 2;
281 ierr = PetscOptionsGetScalarArray(
282 PETSC_NULL,"-param_elastic",paramElastic,&nmax,&flg
283 ); CHKERRQ(ierr);
284
285 paramDensity[0] = 1;
286 paramDensity[1] = 0;
287 int n_max = 2;
288 ierr = PetscOptionsGetScalarArray(PETSC_NULL,"-param_density",paramDensity,&n_max,&flg
289 ); CHKERRQ(ierr);
290
291 ierr = PetscOptionsEnd(); CHKERRQ(ierr);
292
293 PetscFunctionReturn(0);
294}
295
296 /**
297 * returns vertors of data for given index of the point
298 * @param vec_idx reference vector of indices
299 * @param vec_data reference vector of data (colors)
300 * @return error code
301 */
302 PetscErrorCode getDataForGiveIndex(
303 const vector<int> &vec_idx,
304 vector<double> &vec_data
305 ) {
306 PetscFunctionBegin;
307 int size_of_vec_idx = vec_idx.size();
308 vec_data.resize(size_of_vec_idx);
309
310 int ii = 0;
311 for(vector<int>::const_iterator vit = vec_idx.begin();vit!=vec_idx.end();vit++) {
312 vec_data[ii] = metaFile.ElementData(*vit);
313 // if (vec_data[ii] < 500) SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong value (lower than 1000)");
314 // cout << "INDEX "<< *vit << " Data " << vec_data[ii] << endl;
315 ii++;
316 } //cout << "-----------------------" << endl;
317 PetscFunctionReturn(0);
318 }
319
320 /**
321 * returns vertors of coordinates and indices for given coordinate
322 * @param x coordinate of point x
323 * @param y coordinate of point y
324 * @param z coordinate of point z
325 * @param dx dimension of the block in x direction
326 * @param dy dimension of the block in y direction
327 * @param dz dimension of the block in z direction
328 * @param vec_ix reference vector of x coordinates
329 * @param vec_iy reference vector of y coordinates
330 * @param vec_iz reference vector of z coordinates
331 * @param vec_idx reference vector of indices
332 * @return error code
333 */
335 const double x,const double y,const double z, ///< center of the cube
336 const double dx,const double dy,const double dz, ///< size of the cube
337 vector<int> &vec_ix, ///< vectors of coord points inside cube
338 vector<int> &vec_iy, ///<
339 vector<int> &vec_iz, ///<
340 vector<int> &vec_idx, ///< indices of inside cube with dimensions dx dy dz
341 vector<double> &vec_dist ///< vector of squared distances between center of the cube and points
342 ) {
343 PetscFunctionBegin;
344
345 const double *elemet_size = metaFile.ElementSize();
346 const int *dim_size = metaFile.DimSize();
347
348 int nx = ceil(dx/ (elemet_size[0] * sCale)); // number of points within given size of block - dx
349 int ny = ceil(dy/ (elemet_size[1] * sCale)); // number of points within given size of block - dy
350 int nz = ceil(dz/ (elemet_size[2] * sCale)); // number of points within given size of block - dz
351
352 //const double *offset = metaFile.Offset(); // coords before translation (offset)
353 double x_befor_offset = x;
354 double y_befor_offset = y;
355 double z_befor_offset = z;
356 // double x_befor_offset = x - metaFile.Offset()[0]; // probably its not necessary
357 // double y_befor_offset = y - metaFile.Offset()[1];
358 // double z_befor_offset = z - metaFile.Offset()[2];
359 // center of the cube
360 int ix = ceil(x_befor_offset/ (elemet_size[0] * sCale)); // original coords/indices (before multiplying by elem size)
361 int iy = ceil(y_befor_offset/ (elemet_size[1] * sCale));
362 int iz = ceil(z_befor_offset/ (elemet_size[2] * sCale));
363
364 vec_ix.resize(2*nx); // setting sizes of the vectors equal to number of points
365 vec_iy.resize(2*ny);
366 vec_iz.resize(2*nz);
367
368 int ii = 0;
369 for(int i = 0;i<2*nx;i++) { // vectors of coordinates of points within given cube
370 int idx = ix-nx+i;
371 if(idx >= dim_size[0] || idx < 0) continue; // erasing indices beyond the border
372 vec_ix[ii++] = idx;
373 }
374 vec_ix.resize(ii);
375 ii = 0;
376 for(int i = 0;i<2*ny;i++) {
377 int idx = iy-ny+i;
378 if(idx >= dim_size[1] || idx < 0) continue; // erasing indices beyond the border
379 vec_iy[ii++] = idx;
380 }
381 vec_iy.resize(ii);
382 ii = 0;
383 for(int i = 0;i<2*nz;i++) {
384 int idx = iz-nz+i;
385 if(idx >= dim_size[2] || idx < 0) continue; // erasing indices beyond the border
386 vec_iz[ii++] = idx;
387 }
388 vec_iz.resize(ii);
389
390
391 if (dx <= 0) { //
392 vec_ix.resize(1);
393 vec_ix.push_back(round(x_befor_offset/ (elemet_size[0] * sCale)));
394 }
395 if (dy <= 0) {
396 vec_iy.resize(1);
397 vec_iy.push_back(round(y_befor_offset/ (elemet_size[1] * sCale)));
398 }
399 if (dz <= 0) {
400 vec_iz.resize(1);
401 vec_iz.push_back(round(z_befor_offset/ (elemet_size[2] * sCale)));
402 }
403
404 // vec_idx.resize(8*nx*ny*nz);
405 // vec_dist.resize(8*nx*ny*nz);
406 vec_idx.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
407 vec_dist.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
408 int dm_size0_dim_size1 = dim_size[0] * dim_size[1];
409 // calculating the vector of global indices and vector of distances from center of the cube
410
411 if (fIlter == 0) {
412
413 // vector<double> dy2(vec_iy.size());
414 // {
415 // int ii = 0;
416 // for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
417 // double dyy = (*it_iy*elemet_size[1] - y_befor_offset);
418 // dy2[ii++] = dyy * dyy;
419 // }
420 // }
421 // vector<double> dz2(vec_iz.size());
422 // {
423 // int ii = 0;
424 // for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
425 // double dzz = (*it_iz*elemet_size[2] - z_befor_offset);
426 // dz2[ii++] = dzz * dzz;
427 // }
428 // }
429
430 ii = 0;
431 for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
432 double a = (*it_iz) * dm_size0_dim_size1 ;
433 double dz2 = (*it_iz*(elemet_size[2]*sCale) - z_befor_offset) * (*it_iz*(elemet_size[2]*sCale) - z_befor_offset);
434 for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
435 double b = (*it_iy) * dim_size[0];
436 double dy2 = (*it_iy*(elemet_size[1]*sCale) - y_befor_offset) * (*it_iy*(elemet_size[1]*sCale) - y_befor_offset);
437 for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
438 vec_idx[ii] = a + b + (*it_ix);
439 vec_dist[ii] = ( dz2 + dy2 + (*it_ix*(elemet_size[0]*sCale) - x_befor_offset)
440 * (*it_ix*(elemet_size[0]*sCale) - x_befor_offset) );
441 ii++;
442 // cout << "Ix Iy Iz" << vec_idx[ii] << " " << (*it_ix) << " " << (*it_iy) << " " << (*it_iz) << endl;
443 // cout << "Z " << *it_iz << "-" << vec_idx[ii] / (dim_size[0] * dim_size[1]) << " ";
444 // cout << "Y " << *it_iy << "-" << (vec_idx[ii] - (*it_iz) * (dim_size[0] * dim_size[1])) / dim_size[0] << " ";
445 // cout << "X " << *it_ix << "-" << vec_idx[ii] - (*it_iz) * dim_size[0] * dim_size[1] - (*it_iy) * dim_size[0] << " " << endl;
446 // cout << "POINT "<< ii << " GLOBAL " << vec_idx[ii] <<
447 // " DATA "<< metaFile.ElementData(vec_idx[ii]) << " DISTANCE " << vec_dist[ii] << endl;
448
449
450 // if we require more indices than mhd file has ----> ERROR
451 // if(vec_idx[ii] >= metaFile.Quantity()) {
452 // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong index of voxels (> Quantity)");
453 // }
454
455
456 }
457 }
458 }
459
460 }
461 else {
462
463 ii = 0;
464 for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
465 double a = (*it_iz) * dm_size0_dim_size1 ;
466 for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
467 double b = (*it_iy) * dim_size[0];
468 for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
469 vec_idx[ii++] = a + b + (*it_ix);
470 }
471 }
472 }
473 }
474 // cout << ix << " " << " " << iy << " " << iz <<endl;
475 //cout << "--------------------------------------------------------" << endl;
476 PetscFunctionReturn(0);
477 }
478
479 /**
480 * \brief Averages data within given vector
481 *
482 * @param vector data from metaFile
483 * @return averaged data
484 */
485 double getAverage(const vector<double> &vec_data) {
486 PetscFunctionBegin;
487 double data = 0;
488 int i=0;
489 for(vector<double>::const_iterator vec_it = vec_data.begin();vec_it!=vec_data.end();vec_it++) {
490 if (*vec_it < tHreshold) continue;
491 data += *vec_it;
492 i++;
493 }
494 if (i == 0) i = 1; // to avoid division by zero
495 data /= i;
496 //cout <<"data "<< data << " SIZE " << vec_data.size() << endl;
497 return data;
498 }
499
500 /**
501 * \brief Transform data to densities
502 *
503 * @param a array of parameters a,b for linear relation of grayscale - density, a * HU + b = density
504 * @param vec_data data from metaFile
505 * @param vec_density output vector with density
506 * @return vector of densities
507 */
508 PetscErrorCode mapDensityLinear(
509 double *paramDensity,
510 const vector<double> &vec_data,
511 vector<double> &vec_density
512 ) {
513 PetscFunctionBegin;
514 vec_density.resize(vec_data.size());
515 int ii = 0;
516 for(vector<double>::const_iterator vec_it = vec_data.begin();vec_it!=vec_data.end();vec_it++) {
517 vec_density[ii++] = fmax(paramDensity[0] * (*vec_it) + paramDensity[1],0); // // linear density mapping: param_a and param_b
518 //cout << "density " << fmax(a * (*vec_it) + b,0) << endl;
519 //cerr << "tHreshold :" << tHreshold << " fIlter: " << fIlter << " dist toler" << dIst_tol << endl;
520 }
521 PetscFunctionReturn(0);
522 }
523 PetscErrorCode mapDensityLinear(
524 const vector<double> &vec_data,
525 vector<double> &vec_density
526 ) {
527 PetscFunctionBegin;
528 ierr = mapDensityLinear(paramDensity,vec_data,vec_density); CHKERRQ(ierr);
529 PetscFunctionReturn(0);
530 }
531 /**
532 * \brief Transform densities to Youngs modulus
533 *
534 * @param a,b parameter for relation of density - density, b * density ^ a = E
535 * @param rho density
536 * @return youngs modulus
537 */
538 double getElasticMod( // function transforms density data to Young modulus
539 double *param, ///< parameters a,b for density-elasticity relation
540 double rho ///< density
541 ) {
542 return param[1] * pow(rho , param[0]); // ///< power density mapping: param_a and param_b
543 }
544
545 double getElasticMod(double rho) {
546 return getElasticMod(paramElastic,rho); //
547 }
548 /** \brief function returns median of data from given vector
549 *
550 * @param vec_density vector with densities
551 * @return median value
552 */
553 double medianFilter(const vector<double> &vec_density) {
554 vector<double> vec_copy; ///< copy of the vector to allow sorting data
555 vec_copy.resize(vec_density.size());
556 int i = 0;
557 for(vector<double>::const_iterator vec_it = vec_density.begin();vec_it!=vec_density.end();vec_it++){
558 if (*vec_it < tHreshold) continue; ///< applying threshold
559 vec_copy[i++] = *vec_it;
560 //cout << *vec_it << "-";
561 }
562 vec_copy.resize(i);
563
564 //cout << " " << endl;
565 sort(vec_copy.begin(), vec_copy.end()); ///< sorting vector
566 //cout << " " << endl;
567 if (vec_copy.size() % 2 == 0) {
568 return (vec_copy[vec_copy.size() / 2 - 1] + vec_copy[vec_copy.size() / 2]) / 2.0;
569 } else {
570 return vec_copy[vec_copy.size() / 2];
571 }
572 return 0;
573 }
574
575 vector<double> kErnel;
576
577 /**
578 * \brief function smooths values in a given vector depending on distance from center (vec_dist)
579 *
580 */
582 const double sigma, ///< sigma parameter for gausssian smoothing
583 const vector<double> &vec_density, ///< vector with non-smoothed densities
584 const vector<double> &vec_dist ///< vector with squared distances form center of the cubes
585 // vector<double> &vec_density_smooth ///< vector with the smoothed densities
586 ) {
587 //vec_density_smooth.resize(vec_density.size());
588 kErnel.resize(vec_density.size());
589 //calculating kernel
590 int i = 0;
591 double sum = 0;
592 const double k = 2 * sigma*sigma;
593 const double m = (1 / sqrt( M_PI * k));
594 for(int ii = 0; ii < vec_dist.size(); ii++) {
595 if (vec_density[ii] < tHreshold) continue;
596 kErnel[i] = m * exp(- (vec_dist[ii]) / k); // distance is already squared
597 sum += kErnel[i++];
598 //cout <<"KERNEL"<< kernel[i] << endl;;
599 }
600 kErnel.resize(i);
601 ///< normalisation of kernel and smooting values
602 int ii = 0;
603 double smooth_data = 0;
604 for(vector<double>::const_iterator vec_itr = vec_density.begin();vec_itr!=vec_density.end();vec_itr++) {
605 if (*vec_itr < tHreshold) continue;
606 kErnel[ii] /= sum;
607 //cout << kernel[ii] << endl;
608 smooth_data += (*vec_itr) * kErnel[ii++];
609 }
610 if (sum == 0) return 0; ///< avoid division by zero (ad exception) this should not happen /FIXME
611 return smooth_data;
612 }
613
615 const vector<double> &vec_density, ///< vector with non-smoothed densities
616 const vector<double> &vec_dist ///< vector with distances form center of the cubes
617 ) {
618 return mapGaussSmooth(sIgma,vec_density,vec_dist);
619 }
620
621 };
622
623 #endif // WITH_METAIO
624
628 /** \brief Set integrate rule
629 */
630 int getRule(int order) { return 2*order+1; }; // getRule(int order) 2*oRder;
631 };
632
633 /** \brief Calculate volume of the model.
634 * \ingroup bone_remodelling
635 */
637
639 OpVolumeCalculation(const string &field_name,Vec volume_vec):
640 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
641 volumeVec(volume_vec) {
642 }
643
644
645 PetscErrorCode doWork(
646 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
647 ) {
648 PetscFunctionBegin;
649
650 //do it only once, no need to repeat this for edges,faces or tets
651 if(row_type != MBTET) PetscFunctionReturn(0);
652 double vol = getVolume();
653 ierr = VecSetValue(volumeVec,0,vol,ADD_VALUES); CHKERRQ(ierr);
654
655 PetscFunctionReturn(0);
656 }
657
658 };
659 /**
660 * \brief Data structure for storing global and local matrix K and vector f.
661 *
662 \f[
663 Kq=\mathbf{f}
664 \f]
665 *
666 */
667 struct CommonData {
668
669 VectorDouble rhoAtGaussPts; ///< Values of density at integration Pts
670 VectorDouble locF; ///< Local element rhs vector
671 MatrixDouble locA; ///< Local element matrix
672 MatrixDouble transLocA;
673
674 Mat globA; ///< Global matrix
675 Vec globF; ///< Global vector
676
677 boost::shared_ptr<VectorDouble> distanceValue;
678 boost::shared_ptr<MatrixDouble> distanceGradient;
679
680 };
681
682 /** \brief Assemble LHS matrix K.
683 * \ingroup bone_remodelling
684 *
685 \f[
686 K=\int_{V}\phi^{T}\phi\,dV
687 \f]
688 */
690
693
695 const std::string &row_field,const std::string &col_field,CommonData &common_data,DataFromMetaIO &data_from_meta_io
696 ):
697 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
698 row_field,col_field,UserDataOperator::OPROWCOL
699 ),
700 commonData(common_data),
701 dataFromMetaIo(data_from_meta_io) {
702 sYmm = true;
703 }
704
705 PetscErrorCode doWork(
706 int row_side,int col_side,
707 EntityType row_type,EntityType col_type,
708 DataForcesAndSourcesCore::EntData &row_data,
709 DataForcesAndSourcesCore::EntData &col_data
710 ) {
711
712 PetscFunctionBegin;
713
714 int nb_rows = row_data.getIndices().size(); // number of rows
715 int nb_cols = col_data.getIndices().size(); // number of columns
716 if(nb_rows==0) PetscFunctionReturn(0);
717 if(nb_cols==0) PetscFunctionReturn(0);
718
719 commonData.locA.resize(nb_rows,nb_cols,false);
720 commonData.locA.clear();
721 int nb_gauss_pts = row_data.getN().size1();
722 for(int gg = 0;gg<nb_gauss_pts;gg++) {
723 double vol = getVolume()*getGaussPts()(3,gg);
724 VectorAdaptor row_base_function = row_data.getN(gg);
725 VectorAdaptor col_base_function = col_data.getN(gg);
726 commonData.locA += vol*outer_prod(row_base_function,col_base_function);
727 MatrixAdaptor row_diff_base_function = row_data.getDiffN(gg);
728 MatrixAdaptor col_diff_base_function = col_data.getDiffN(gg);
729
730 //const double lambda = 1;
732 row_diff_base_function,trans(col_diff_base_function));
733 }
734 ierr = MatSetValues(
736 row_data.getIndices().size(),
737 &row_data.getIndices()[0],
738 col_data.getIndices().size(),
739 &col_data.getIndices()[0],
740 &commonData.locA(0,0),
741 ADD_VALUES
742 ); CHKERRQ(ierr);
743
744 // Assemble off symmetric side
745 if(row_side!=col_side || row_type!=col_type) {
746 commonData.transLocA.resize(commonData.locA.size2(),commonData.locA.size1(),false);
747 noalias(commonData.transLocA) = trans(commonData.locA);
748 ierr = MatSetValues(
750 col_data.getIndices().size(),
751 &col_data.getIndices()[0],
752 row_data.getIndices().size(),
753 &row_data.getIndices()[0],
754 &commonData.transLocA(0,0),
755 ADD_VALUES
756 ); CHKERRQ(ierr);
757 }
758
759 PetscFunctionReturn(0);
760 }
761
762 };
763
764 #ifdef WITH_METAIO
765
766 /** \brief Assemble local vector containing density data.
767 * \ingroup bone_remodelling
768 */
770
774
776 const string &row_field,
777 CommonData &common_data,
778 DataFromMetaIO &data_from_meta_io,
779 SurfaceKDTree &surface_distance
780 ):
781 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
782 row_field,row_field,UserDataOperator::OPROW
783 ),
784 commonData(common_data),
785 dataFromMetaIo(data_from_meta_io),
786 surfaceDistance(surface_distance) {
787 }
788
789
790 PetscErrorCode doWork(
791 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
792 ) {
793 PetscFunctionBegin;
794
795 if(row_type!=MBVERTEX) PetscFunctionReturn(0);
796
797 int nb_rows = row_data.getIndices().size();
798 if(nb_rows == 0) PetscFunctionReturn(0);
799
800 ////////////////////////////////////////////////////////////////////////////////////////////
801 vector<int> vecIx; ///< reference vectors of coordinates
802 vector<int> vecIy;
803 vector<int> vecIz;
804 vector<int> vecIdx; ///< reference vector of indices
805
806 vector<double> vecData; ///< reference vector of data (colors)
807 vector<double> vecDensity; ///< vector of the densities
808 vector<double> vecDist; ///< vector of reference distances of the points
809 ///////////////////////////////////////////////////////////////////////////////////////////////////
810
811 auto t0_dist = getFTensor0FromVec(*commonData.distanceValue);
812 auto t1_grad = getFTensor1FromMat<3>(*commonData.distanceGradient);
814
816 &getCoordsAtGaussPts()(0,0),
817 &getCoordsAtGaussPts()(0,1),
818 &getCoordsAtGaussPts()(0,2),
819 3
820 );
821
823 dataFromMetaIo.metaFile.ElementSize()[0] * dataFromMetaIo.sCale,
824 dataFromMetaIo.metaFile.ElementSize()[1] * dataFromMetaIo.sCale,
826 );
827 t1_dx(i) *= dataFromMetaIo.cUbe_size;
828 double cube_size = sqrt(t1_dx(i)*t1_dx(i))*dataFromMetaIo.dIst_tol;
829
830 const int nb_gauss_pts = row_data.getN().size1();
831 commonData.rhoAtGaussPts.resize(nb_gauss_pts,false);
832 for(int gg = 0;gg<nb_gauss_pts;gg++) {
833
834 if(t0_dist<cube_size) {
835 const double nrm2_grad = sqrt(t1_grad(i)*t1_grad(i));
837 t1_delta(i) = t1_grad(i)/nrm2_grad;
838 t1_xyz(i) -= (t0_dist-cube_size)*t1_delta(i);
839 }
840
841
842 // //double cube_size = sqrt(dx*dx+dy*dy+dz*dz)*dataFromMetaIo.dIst_tol; //FIXME
843 // double xc_dm = (2 * dataFromMetaIo.cUbe_size - 1) * dataFromMetaIo.metaFile.ElementSize()[0];
844 // double yc_dm = (2 * dataFromMetaIo.cUbe_size - 1) * dataFromMetaIo.metaFile.ElementSize()[1];
845 // double zc_dm = (2 * dataFromMetaIo.cUbe_size - 1) * dataFromMetaIo.metaFile.ElementSize()[2];
846 // double cube_size = 0.5 * sqrt(xc_dm * xc_dm + yc_dm * yc_dm + zc_dm * zc_dm) * dataFromMetaIo.dIst_tol; // real cube size /2
847 //
848 // if(dist<cube_size) {
849 // // Shift cube from boundary
850 // x += (dist-cube_size)*delta_x/dist;
851 // y += (dist-cube_size)*delta_y/dist;
852 // z += (dist-cube_size)*delta_z/dist;
853 // }
854
855 //Calculate density like in OpMassCalculation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
857 t1_xyz(0),t1_xyz(1),t1_xyz(2),
858 t1_dx(0),t1_dx(1),t1_dx(2),
859 vecIx,vecIy,vecIz,vecIdx, vecDist
860 ); CHKERRQ(ierr);
861
862 ierr = dataFromMetaIo.getDataForGiveIndex(vecIdx,vecData); CHKERRQ(ierr);
863 ierr = dataFromMetaIo.mapDensityLinear(vecData,vecDensity); CHKERRQ(ierr);
864
865 double data = 0.0;
866 switch (dataFromMetaIo.fIlter) {
867 case 0:
868 data = dataFromMetaIo.mapGaussSmooth(vecDensity,vecDist);
869 break;
870 case 1:
871 data = dataFromMetaIo.getAverage(vecDensity);
872 break;
873 case 2:
874 data = dataFromMetaIo.medianFilter(vecDensity);
875
876 default:
877 data = dataFromMetaIo.getAverage(vecDensity);
878 }
879
880 commonData.rhoAtGaussPts[gg] = data;//pow(x,2)*y; //data; // not multiply this by vol
881
882 ++t1_xyz;
883 ++t0_dist;
884 ++t1_grad;
885
886 }
887
888 PetscFunctionReturn(0);
889 }
890
891
892 };
893
894 /**
895 * \brief Assemble RHS vector f.
896 * \ingroup bone_remodelling
897 *
898 \f[
899 \mathbf{f}_{e}=\int_{V_{e}}\phi q^{\ast}\,dV_{e}
900 \f]
901 */
903
905
907 const string &row_field,
908 CommonData &common_data
909 ):
910 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
911 row_field,row_field,UserDataOperator::OPROW
912 ),
913 commonData(common_data) {
914 }
915
916
917 PetscErrorCode doWork(
918 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
919 ) {
920 PetscFunctionBegin;
921
922 int nb_rows = row_data.getIndices().size();
923 if(nb_rows == 0) PetscFunctionReturn(0);
924 commonData.locF.resize(nb_rows,false);
925 commonData.locF.clear();
926
927 int nb_gauss_pts = row_data.getN().size1();
928 for(int gg = 0;gg<nb_gauss_pts;gg++) {
929
930 double vol = getVolume()*getGaussPts()(3,gg);
931 VectorAdaptor base_function = row_data.getN(gg);
932 commonData.locF += vol*commonData.rhoAtGaussPts[gg]*base_function;
933
934 }
935
936 ierr = VecSetValues(
938 row_data.getIndices().size(),
939 &row_data.getIndices()[0],
940 &commonData.locF[0],
941 ADD_VALUES
942 ); CHKERRQ(ierr);
943
944 PetscFunctionReturn(0);
945 }
946
947 };
948
949 /**
950 * \brief Calculate mass before approximation.
951 * \ingroup bone_remodelling
952 *
953 \f[
954 M=\int_{V}\rho\,dV
955 \f]
956 */
958
961
963 const string &field_name,Vec mass_vec,DataFromMetaIO &data_from_meta_io
964 ):
965 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
966 massVec(mass_vec),
967 dataFromMetaIo(data_from_meta_io) {
968 }
969
970 vector<int> vecIx; ///< reference vectors of coordinates
971 vector<int> vecIy;
972 vector<int> vecIz;
973 vector<int> vecIdx; ///< reference vector of indices
974
975 vector<double> vecData; ///< reference vector of data (colors)
976 vector<double> vecDensity; ///< vector of the densities
977 vector<double> vecDist; ///< vector of reference distances of the points
978
979
980 PetscErrorCode doWork(
981 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
982 ) {
983 PetscFunctionBegin;
984
985 //do it only once, no need to repeat this for edges,faces or tets
986 if(row_type != MBVERTEX) PetscFunctionReturn(0);
987 //double dx, dy, dz; dx = dy = dz = dataFromMetaIo.cUbe_size * dataFromMetaIo.metaFile.ElementSize()[0];
988
989 double dx = dataFromMetaIo.cUbe_size * dataFromMetaIo.metaFile.ElementSize()[0] * dataFromMetaIo.sCale;
990 double dy = dataFromMetaIo.cUbe_size * dataFromMetaIo.metaFile.ElementSize()[1] * dataFromMetaIo.sCale;
991 double dz = dataFromMetaIo.cUbe_size * dataFromMetaIo.metaFile.ElementSize()[2] * dataFromMetaIo.sCale;
992 //cout << " DX DY DZ " << dx << " " << dy << " " << dz << endl;
993 int nb_gauss_pts = row_data.getN().size1();
994 for(int gg = 0;gg<nb_gauss_pts;gg++) {
995
996 double vol = getVolume()*getGaussPts()(3,gg);
997 double x = getCoordsAtGaussPts()(gg,0);
998 double y = getCoordsAtGaussPts()(gg,1);
999 double z = getCoordsAtGaussPts()(gg,2);
1000 //cout <<"i: " << gg << " x " << x << " y " << y << " z " << z <<endl;
1001
1003 x,y,z,
1004 dx,dy,dz,
1006 ); CHKERRQ(ierr);
1007
1010
1011 //double data = dataFromMetaIo.mapGaussSmooth(vecDensity,vecDist);
1012 //double data = dataFromMetaIo.medianFilter(vecDensity);
1013 //double data = dataFromMetaIo.getAverage(vecDensity);
1014
1015 double data = 0.0;
1016 switch (dataFromMetaIo.fIlter)
1017 {
1018 case 0:
1020 break;
1021 case 1:
1023 break;
1024 case 2:
1026
1027 default:
1029 }
1030
1031
1032 double rho = data * vol;
1033 ierr = VecSetValue(massVec,0,rho,ADD_VALUES); CHKERRQ(ierr);
1034
1035 }
1036
1037 PetscFunctionReturn(0);
1038 }
1039
1040 };
1041
1042 /**
1043 * \brief Calculate mass after approximation.
1044 * \ingroup bone_remodelling
1045 *
1046 \f[
1047 M=\int_{V}\rho N\,dV
1048 \f]
1049 */
1051
1053
1054
1056 const string &field_name,Vec mass_vec_approx
1057 ):
1058 MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
1059 massVec_approx(mass_vec_approx){
1060 }
1061
1062
1063 PetscErrorCode doWork(
1064 int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1065 ) {
1066 PetscFunctionBegin;
1067
1068 int nb_rows = row_data.getIndices().size();
1069 if(!nb_rows) {
1070 PetscFunctionReturn(0);
1071 }
1072
1073 int nb_gauss_pts = row_data.getN().size1();
1074
1075 for(int gg = 0;gg<nb_gauss_pts;gg++) {
1076
1077 double vol = getVolume()*getGaussPts()(3,gg);
1078
1079 VectorAdaptor shape_function = row_data.getN(gg,nb_rows);
1080 double rho = 0;
1081 rho = inner_prod(shape_function,row_data.getFieldData());
1082 //cout << rho << " " << shape_function << " " << row_data.getFieldData() << endl;
1083
1084 double mass = rho * vol;
1085 ierr = VecSetValue(massVec_approx,0,mass,ADD_VALUES); CHKERRQ(ierr);
1086
1087 }
1088 //cout << endl;
1089 PetscFunctionReturn(0);
1090 }
1091
1092 };
1093 #endif // WITH_METAIO
1094
1095
1096}
1097
1098#endif //__DENSITY_MAPS_HPP__
1099
1100/***************************************************************************//**
1101* \defgroup bone_remodelling Bone remodeling
1102* \ingroup user_modules
1103*
1104* \brief Bone remodeling module
1105*
1106******************************************************************************/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
static PetscErrorCode ierr
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
constexpr int order
FTensor::Index< 'm', SPACE_DIM > m
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'k', 3 > k
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double rho
Definition: plastic.cpp:187
constexpr auto field_name
static constexpr double delta
Data structure for storing global and local matrix K and vector f.
Vec globF
Global vector.
VectorDouble rhoAtGaussPts
Values of density at integration Pts.
MatrixDouble locA
Local element matrix.
boost::shared_ptr< MatrixDouble > distanceGradient
Mat globA
Global matrix.
VectorDouble locF
Local element rhs vector.
boost::shared_ptr< VectorDouble > distanceValue
Load data from MetaImage file, translate grayscale values to densities.
double mapGaussSmooth(const vector< double > &vec_density, const vector< double > &vec_dist)
PetscErrorCode getDataForGiveIndex(const vector< int > &vec_idx, vector< double > &vec_data)
double medianFilter(const vector< double > &vec_density)
function returns median of data from given vector
DataFromMetaIO(MetaImage &metaFile)
double getElasticMod(double rho)
PetscErrorCode mapDensityLinear(double *paramDensity, const vector< double > &vec_data, vector< double > &vec_density)
Transform data to densities.
PetscErrorCode getInidcesForGivenCoordinate(const double x, const double y, const double z, const double dx, const double dy, const double dz, vector< int > &vec_ix, vector< int > &vec_iy, vector< int > &vec_iz, vector< int > &vec_idx, vector< double > &vec_dist)
PetscErrorCode mapDensityLinear(const vector< double > &vec_data, vector< double > &vec_density)
double getAverage(const vector< double > &vec_data)
Averages data within given vector.
double getElasticMod(double *param, double rho)
Transform densities to Youngs modulus.
double mapGaussSmooth(const double sigma, const vector< double > &vec_density, const vector< double > &vec_dist)
function smooths values in a given vector depending on distance from center (vec_dist)
DensityMapFe(MoFEM::Interface &m_field)
int getRule(int order)
Set integrate rule.
Assemble RHS vector f.
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpAssmbleRhs(const string &row_field, CommonData &common_data)
Assemble LHS matrix K.
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpCalculateLhs(const std::string &row_field, const std::string &col_field, CommonData &common_data, DataFromMetaIO &data_from_meta_io)
Assemble local vector containing density data.
OpCalulatefRhoAtGaussPts(const string &row_field, CommonData &common_data, DataFromMetaIO &data_from_meta_io, SurfaceKDTree &surface_distance)
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Calculate mass after approximation.
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpMassCalculationFromApprox(const string &field_name, Vec mass_vec_approx)
Calculate mass before approximation.
OpMassCalculation(const string &field_name, Vec mass_vec, DataFromMetaIO &data_from_meta_io)
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
vector< int > vecIx
reference vectors of coordinates
vector< double > vecDensity
vector of the densities
vector< double > vecDist
vector of reference distances of the points
vector< int > vecIdx
reference vector of indices
vector< double > vecData
reference vector of data (colors)
Calculate volume of the model.
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
OpVolumeCalculation(const string &field_name, Vec volume_vec)
Create KDTree on surface of the mesh and calculate distance.
Definition: DensityMaps.hpp:32
PetscErrorCode findClosestPointToTheSurface(const double x, const double y, const double z, double &dx, double &dy, double &dz)
Find point on surface which is closet to given point coordinates.
Definition: DensityMaps.hpp:88
SurfaceKDTree(MoFEM::Interface &m_field)
Definition: DensityMaps.hpp:39
PetscErrorCode setDistanceFromSurface(const std::string field_name)
Set distance to vertices on mesh.
PetscErrorCode takeASkin(Range &volume)
Take a skin from a mesh.
Definition: DensityMaps.hpp:55
PetscErrorCode buildTree()
Build tree.
Definition: DensityMaps.hpp:72
MoFEM::Interface & mField
Definition: DensityMaps.hpp:34
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)