v0.14.0
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 
27 namespace BoneRemodeling {
28  /**
29  * \brief Create KDTree on surface of the mesh and calculate distance.
30  * \ingroup bone_remodelling
31  */
32  struct SurfaceKDTree {
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;
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  */
205  struct DataFromMetaIO {
206 
207  MetaImage &metaFile;
209  metaFile(metaFile) {
210  }
211 
212 
213  double sIgma;
214  double cUbe_size;
215  int fIlter;
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;
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
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,
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 
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,
825  dataFromMetaIo.metaFile.ElementSize()[2] * 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 ******************************************************************************/
BoneRemodeling::OpVolumeCalculation::OpVolumeCalculation
OpVolumeCalculation(const string &field_name, Vec volume_vec)
Definition: DensityMaps.hpp:639
BoneRemodeling::OpVolumeCalculation::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DensityMaps.hpp:645
BoneRemodeling::OpMassCalculation::vecIdx
vector< int > vecIdx
reference vector of indices
Definition: DensityMaps.hpp:973
BoneRemodeling::OpVolumeCalculation::volumeVec
Vec volumeVec
Definition: DensityMaps.hpp:638
BoneRemodeling::DataFromMetaIO::DataFromMetaIO
DataFromMetaIO(MetaImage &metaFile)
Definition: DensityMaps.hpp:208
BoneRemodeling::OpMassCalculation::vecIx
vector< int > vecIx
reference vectors of coordinates
Definition: DensityMaps.hpp:970
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
BoneRemodeling::OpCalulatefRhoAtGaussPts::commonData
CommonData & commonData
Definition: DensityMaps.hpp:771
MoFEM::PetscOptionsGetScalarArray
PetscErrorCode PetscOptionsGetScalarArray(PetscOptions options, const char pre[], const char name[], PetscScalar dvalue[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:228
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
rho
double rho
Definition: plastic.cpp:140
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
BoneRemodeling::OpMassCalculationFromApprox::OpMassCalculationFromApprox
OpMassCalculationFromApprox(const string &field_name, Vec mass_vec_approx)
Definition: DensityMaps.hpp:1055
BoneRemodeling::CommonData::locF
VectorDouble locF
Local element rhs vector.
Definition: DensityMaps.hpp:670
BoneRemodeling::OpMassCalculation::vecDist
vector< double > vecDist
vector of reference distances of the points
Definition: DensityMaps.hpp:977
BoneRemodeling::OpMassCalculation::vecIz
vector< int > vecIz
Definition: DensityMaps.hpp:972
BoneRemodeling::OpMassCalculationFromApprox::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DensityMaps.hpp:1063
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
BoneRemodeling::CommonData::rhoAtGaussPts
VectorDouble rhoAtGaussPts
Values of density at integration Pts.
Definition: DensityMaps.hpp:669
BoneRemodeling::DataFromMetaIO::getElasticMod
double getElasticMod(double rho)
Definition: DensityMaps.hpp:545
BoneRemodeling::OpMassCalculationFromApprox::massVec_approx
Vec massVec_approx
Definition: DensityMaps.hpp:1052
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
BoneRemodeling::DataFromMetaIO::cUbe_size
double cUbe_size
Definition: DensityMaps.hpp:214
BoneRemodeling::DataFromMetaIO::paramDensity
double paramDensity[2]
Definition: DensityMaps.hpp:220
BoneRemodeling::OpMassCalculation::OpMassCalculation
OpMassCalculation(const string &field_name, Vec mass_vec, DataFromMetaIO &data_from_meta_io)
Definition: DensityMaps.hpp:962
BoneRemodeling::OpCalulatefRhoAtGaussPts
Assemble local vector containing density data.
Definition: DensityMaps.hpp:769
BoneRemodeling::OpMassCalculationFromApprox
Calculate mass after approximation.
Definition: DensityMaps.hpp:1050
BoneRemodeling::OpCalulatefRhoAtGaussPts::surfaceDistance
SurfaceKDTree & surfaceDistance
Definition: DensityMaps.hpp:773
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
BoneRemodeling::SurfaceKDTree::kdTreeRootMeshset
EntityHandle kdTreeRootMeshset
Definition: DensityMaps.hpp:35
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
BoneRemodeling::OpCalculateLhs::OpCalculateLhs
OpCalculateLhs(const std::string &row_field, const std::string &col_field, CommonData &common_data, DataFromMetaIO &data_from_meta_io)
Definition: DensityMaps.hpp:694
BoneRemodeling::DataFromMetaIO::mapDensityLinear
PetscErrorCode mapDensityLinear(double *paramDensity, const vector< double > &vec_data, vector< double > &vec_density)
Transform data to densities.
Definition: DensityMaps.hpp:508
order
constexpr int order
Definition: dg_projection.cpp:18
BoneRemodeling::OpMassCalculation::massVec
Vec massVec
Definition: DensityMaps.hpp:959
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
BoneRemodeling::CommonData
Data structure for storing global and local matrix K and vector f.
Definition: DensityMaps.hpp:667
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1235
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
BoneRemodeling::DensityMapFe::DensityMapFe
DensityMapFe(MoFEM::Interface &m_field)
Definition: DensityMaps.hpp:626
BoneRemodeling::DataFromMetaIO::getElasticMod
double getElasticMod(double *param, double rho)
Transform densities to Youngs modulus.
Definition: DensityMaps.hpp:538
BoneRemodeling::SurfaceKDTree::SurfaceKDTree
SurfaceKDTree(MoFEM::Interface &m_field)
Definition: DensityMaps.hpp:39
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
BoneRemodeling::OpAssmbleRhs::commonData
CommonData & commonData
Definition: DensityMaps.hpp:904
BoneRemodeling::SurfaceKDTree
Create KDTree on surface of the mesh and calculate distance.
Definition: DensityMaps.hpp:32
BoneRemodeling::SurfaceKDTree::buildTree
PetscErrorCode buildTree()
Build tree.
Definition: DensityMaps.hpp:72
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
BoneRemodeling::OpCalculateLhs::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: DensityMaps.hpp:705
BoneRemodeling::DataFromMetaIO::sCale
double sCale
Definition: DensityMaps.hpp:216
BoneRemodeling::OpVolumeCalculation
Calculate volume of the model.
Definition: DensityMaps.hpp:636
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
BoneRemodeling::DataFromMetaIO::fromOptions
PetscErrorCode fromOptions()
Definition: DensityMaps.hpp:223
a
constexpr double a
Definition: approx_sphere.cpp:30
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
BoneRemodeling::DataFromMetaIO::lAmbda
double lAmbda
Definition: DensityMaps.hpp:221
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
BoneRemodeling::OpAssmbleRhs
Assemble RHS vector f.
Definition: DensityMaps.hpp:902
BoneRemodeling::CommonData::distanceValue
boost::shared_ptr< VectorDouble > distanceValue
Definition: DensityMaps.hpp:677
convert.type
type
Definition: convert.py:64
BoneRemodeling::DensityMapFe::getRule
int getRule(int order)
Set integrate rule.
Definition: DensityMaps.hpp:630
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
BoneRemodeling::OpCalulatefRhoAtGaussPts::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DensityMaps.hpp:790
BoneRemodeling::DataFromMetaIO::kErnel
vector< double > kErnel
Definition: DensityMaps.hpp:575
BoneRemodeling::SurfaceKDTree::kdTree
AdaptiveKDTree kdTree
Definition: DensityMaps.hpp:36
BoneRemodeling
Definition: DensityMaps.hpp:27
BoneRemodeling::OpCalculateLhs::dataFromMetaIo
DataFromMetaIO & dataFromMetaIo
Definition: DensityMaps.hpp:692
BoneRemodeling::DataFromMetaIO::getInidcesForGivenCoordinate
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)
Definition: DensityMaps.hpp:334
BoneRemodeling::OpCalulatefRhoAtGaussPts::dataFromMetaIo
DataFromMetaIO & dataFromMetaIo
Definition: DensityMaps.hpp:772
BoneRemodeling::DataFromMetaIO::getDataForGiveIndex
PetscErrorCode getDataForGiveIndex(const vector< int > &vec_idx, vector< double > &vec_data)
Definition: DensityMaps.hpp:302
MoFEM::CoreInterface::get_field_ents
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
BoneRemodeling::SurfaceKDTree::mField
MoFEM::Interface & mField
Definition: DensityMaps.hpp:34
MoFEM::ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
Definition: ForcesAndSourcesCore.hpp:1264
BoneRemodeling::SurfaceKDTree::sKin
Range sKin
Definition: DensityMaps.hpp:47
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
BoneRemodeling::OpCalulatefRhoAtGaussPts::OpCalulatefRhoAtGaussPts
OpCalulatefRhoAtGaussPts(const string &row_field, CommonData &common_data, DataFromMetaIO &data_from_meta_io, SurfaceKDTree &surface_distance)
Definition: DensityMaps.hpp:775
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
BoneRemodeling::DataFromMetaIO::dIst_tol
double dIst_tol
Definition: DensityMaps.hpp:217
BoneRemodeling::DataFromMetaIO::medianFilter
double medianFilter(const vector< double > &vec_density)
function returns median of data from given vector
Definition: DensityMaps.hpp:553
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
BoneRemodeling::CommonData::distanceGradient
boost::shared_ptr< MatrixDouble > distanceGradient
Definition: DensityMaps.hpp:678
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
BoneRemodeling::DataFromMetaIO::sIgma
double sIgma
Definition: DensityMaps.hpp:213
BoneRemodeling::DataFromMetaIO::fIlter
int fIlter
Definition: DensityMaps.hpp:215
BoneRemodeling::DataFromMetaIO::mapDensityLinear
PetscErrorCode mapDensityLinear(const vector< double > &vec_data, vector< double > &vec_density)
Definition: DensityMaps.hpp:523
BoneRemodeling::CommonData::globF
Vec globF
Global vector.
Definition: DensityMaps.hpp:675
BoneRemodeling::SurfaceKDTree::takeASkin
PetscErrorCode takeASkin(Range &volume)
Take a skin from a mesh.
Definition: DensityMaps.hpp:55
BoneRemodeling::OpMassCalculation::vecDensity
vector< double > vecDensity
vector of the densities
Definition: DensityMaps.hpp:976
Range
BoneRemodeling::OpCalculateLhs::commonData
CommonData & commonData
Definition: DensityMaps.hpp:691
LOBATTO_PHI0
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
Definition: base_functions.h:126
BoneRemodeling::OpCalculateLhs
Assemble LHS matrix K.
Definition: DensityMaps.hpp:689
BoneRemodeling::OpMassCalculation::vecData
vector< double > vecData
reference vector of data (colors)
Definition: DensityMaps.hpp:975
BoneRemodeling::DataFromMetaIO
Load data from MetaImage file, translate grayscale values to densities.
Definition: DensityMaps.hpp:205
BoneRemodeling::CommonData::transLocA
MatrixDouble transLocA
Definition: DensityMaps.hpp:672
BoneRemodeling::OpMassCalculation::vecIy
vector< int > vecIy
Definition: DensityMaps.hpp:971
BoneRemodeling::DataFromMetaIO::paramElastic
double paramElastic[2]
Definition: DensityMaps.hpp:219
BoneRemodeling::OpMassCalculation
Calculate mass before approximation.
Definition: DensityMaps.hpp:957
BoneRemodeling::OpMassCalculation::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DensityMaps.hpp:980
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
BoneRemodeling::SurfaceKDTree::setDistanceFromSurface
PetscErrorCode setDistanceFromSurface(const std::string field_name)
Set distance to vertices on mesh.
Definition: DensityMaps.hpp:119
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
BoneRemodeling::DataFromMetaIO::mapGaussSmooth
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)
Definition: DensityMaps.hpp:581
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
BoneRemodeling::CommonData::locA
MatrixDouble locA
Local element matrix.
Definition: DensityMaps.hpp:671
BoneRemodeling::OpAssmbleRhs::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: DensityMaps.hpp:917
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::Types::MatrixAdaptor
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:132
BoneRemodeling::CommonData::globA
Mat globA
Global matrix.
Definition: DensityMaps.hpp:674
BoneRemodeling::OpMassCalculation::dataFromMetaIo
DataFromMetaIO & dataFromMetaIo
Definition: DensityMaps.hpp:960
BoneRemodeling::DataFromMetaIO::mapGaussSmooth
double mapGaussSmooth(const vector< double > &vec_density, const vector< double > &vec_dist)
Definition: DensityMaps.hpp:614
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
BoneRemodeling::DataFromMetaIO::getAverage
double getAverage(const vector< double > &vec_data)
Averages data within given vector.
Definition: DensityMaps.hpp:485
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
BoneRemodeling::OpAssmbleRhs::OpAssmbleRhs
OpAssmbleRhs(const string &row_field, CommonData &common_data)
Definition: DensityMaps.hpp:906
BoneRemodeling::SurfaceKDTree::findClosestPointToTheSurface
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
BoneRemodeling::DataFromMetaIO::tHreshold
double tHreshold
Definition: DensityMaps.hpp:218
BoneRemodeling::DensityMapFe
Definition: DensityMaps.hpp:625
BoneRemodeling::DataFromMetaIO::metaFile
MetaImage & metaFile
Definition: DensityMaps.hpp:207
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567