v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
BoneRemodeling::SurfaceKDTree Struct Reference
User modules » Bone remodeling

Create KDTree on surface of the mesh and calculate distance. More...

#include <users_modules/bone_remodelling/src/DensityMaps.hpp>

Collaboration diagram for BoneRemodeling::SurfaceKDTree:
[legend]

Public Member Functions

 SurfaceKDTree (MoFEM::Interface &m_field)
 
PetscErrorCode takeASkin (Range &volume)
 Take a skin from a mesh. More...
 
PetscErrorCode buildTree ()
 Build tree. More...
 
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. More...
 
PetscErrorCode setDistanceFromSurface (const std::string field_name)
 Set distance to vertices on mesh. More...
 

Public Attributes

MoFEM::Interface & mField
 
EntityHandle kdTreeRootMeshset
 
AdaptiveKDTree kdTree
 
Range sKin
 

Detailed Description

Create KDTree on surface of the mesh and calculate distance.

Definition at line 32 of file DensityMaps.hpp.

Constructor & Destructor Documentation

◆ SurfaceKDTree()

BoneRemodeling::SurfaceKDTree::SurfaceKDTree ( MoFEM::Interface &  m_field)
inline

Definition at line 39 of file DensityMaps.hpp.

39 :
40 mField(m_field),
41 kdTree(&m_field.get_moab()) {
42 }
BoneRemodeling::SurfaceKDTree::mField
MoFEM::Interface & mField
Definition: DensityMaps.hpp:34
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0

Member Function Documentation

◆ buildTree()

PetscErrorCode BoneRemodeling::SurfaceKDTree::buildTree ( )
inline

Build tree.

Returns
Error code

Definition at line 72 of file DensityMaps.hpp.

72 {
73 PetscFunctionBegin;
74 rval = kdTree.build_tree(sKin,&kdTreeRootMeshset); CHKERRQ_MOAB(rval);
75 PetscFunctionReturn(0);
76 }
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454

◆ findClosestPointToTheSurface()

PetscErrorCode BoneRemodeling::SurfaceKDTree::findClosestPointToTheSurface ( const double  x,
const double  y,
const double  z,
double &  dx,
double &  dy,
double &  dz 
)
inline

Find point on surface which is closet to given point coordinates.

Parameters
x1st coordinate of point
y2nd coordinate of point
z3rd coordinate of point
dxcoordinate of distance vector
dycoordinate of distance vector
zcoordinate of distance vector
Returns
Error code

Definition at line 88 of file DensityMaps.hpp.

95 {
96 PetscFunctionBegin;
97 const double coords[] = {x,y,z};
98 double closest_point[3];
99 EntityHandle triangle;
100 rval = kdTree.closest_triangle(
101 kdTreeRootMeshset,
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 }

◆ setDistanceFromSurface()

PetscErrorCode BoneRemodeling::SurfaceKDTree::setDistanceFromSurface ( const std::string  field_name)
inline

Set distance to vertices on mesh.

Parameters
field_nameName of field for distances
Returns
Error Code

Definition at line 119 of file DensityMaps.hpp.

119 {
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];
137 ierr = findClosestPointToTheSurface(
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];
162 ierr = findClosestPointToTheSurface(
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 }
ierr
static PetscErrorCode ierr
Definition: base_functions.c:13
LOBATTO_PHI0
#define LOBATTO_PHI0(x)
Definitions taken from Hermes2d code.
Definition: base_functions.h:126
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::CoreInterface::get_field_ents
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
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
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number

◆ takeASkin()

PetscErrorCode BoneRemodeling::SurfaceKDTree::takeASkin ( Range &  volume)
inline

Take a skin from a mesh.

Parameters
volumeVolume elements in the mesh
Returns
Error code

Definition at line 55 of file DensityMaps.hpp.

55 {
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 }

Member Data Documentation

◆ kdTree

AdaptiveKDTree BoneRemodeling::SurfaceKDTree::kdTree

Definition at line 36 of file DensityMaps.hpp.

◆ kdTreeRootMeshset

EntityHandle BoneRemodeling::SurfaceKDTree::kdTreeRootMeshset

Definition at line 35 of file DensityMaps.hpp.

◆ mField

MoFEM::Interface& BoneRemodeling::SurfaceKDTree::mField

Definition at line 34 of file DensityMaps.hpp.

◆ sKin

Range BoneRemodeling::SurfaceKDTree::sKin

Definition at line 47 of file DensityMaps.hpp.

The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.5 and hosted at University of Glasgow