11 #if !defined(KRATOS_TETGEN_PFEM_MODELER_VMS_H_INCLUDED )
12 #define KRATOS_TETGEN_PFEM_MODELER_VMS_H_INCLUDED
20 #include <boost/timer.hpp>
28 #include "utilities/geometry_utilities.h"
105 Element const& rReferenceElement,
106 Condition const& rReferenceBoundaryCondition,
108 double alpha_param = 1.4,
double h_factor=0.5)
112 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_FREE_SURFACE)==
false )
113 KRATOS_THROW_ERROR(std::logic_error,
"Add ----IS_FREE_SURFACE---- variable!!!!!! ERROR",
"");
114 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_STRUCTURE)==
false )
115 KRATOS_THROW_ERROR(std::logic_error,
"Add ----IS_STRUCTURE---- variable!!!!!! ERROR",
"");
116 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_BOUNDARY)==
false )
117 KRATOS_THROW_ERROR(std::logic_error,
"Add ----IS_BOUNDARY---- variable!!!!!! ERROR",
"");
118 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_FLUID)==
false )
121 KRATOS_WATCH(
" ENTERED TETGENMESHSUITE PFEM of Meshing Application")
128 boost::timer auxiliary;
131 typedef Node::Pointer PointPointerType;
133 typedef std::vector<PointType::Pointer>
PointVector;
151 unsigned int bucket_size = 20;
154 unsigned int max_results = 50;
159 Node work_point(0,0.0,0.0,0.0);
166 list_of_nodes.reserve(ThisModelPart.
Nodes().size());
167 for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.
NodesBegin() ; i_node != ThisModelPart.
NodesEnd() ; i_node++)
169 (list_of_nodes).push_back(*(i_node.base()));
173 kd_tree nodes_tree1(list_of_nodes.begin(),list_of_nodes.end(), bucket_size);
176 unsigned int n_points_in_radius;
180 for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.
NodesBegin(); in != ThisModelPart.
NodesEnd(); in++)
182 radius=h_factor*in->FastGetSolutionStepValue(NODAL_H);
184 work_point[0]=in->
X();
185 work_point[1]=in->
Y();
186 work_point[2]=in->
Z();
188 n_points_in_radius = nodes_tree1.SearchInRadius(work_point,
radius,
res.begin(),res_distances.begin(), max_results);
189 if (n_points_in_radius>1)
191 if (in->FastGetSolutionStepValue(IS_BOUNDARY)==0.0 && in->FastGetSolutionStepValue(IS_STRUCTURE)==0.0)
194 double erased_nodes = 0;
195 for(
auto i=
res.begin();
i!=
res.begin() + n_points_in_radius ;
i++)
196 erased_nodes += in->Is(TO_ERASE);
198 if( erased_nodes < 1)
199 in->Set(TO_ERASE,
true);
201 else if ( (in)->FastGetSolutionStepValue(IS_STRUCTURE)!=1.0)
207 for(
auto i=
res.begin();
i!=
res.begin() + n_points_in_radius ;
i++)
209 if ( (*i)->FastGetSolutionStepValue(IS_BOUNDARY,1)==1.0 && res_distances[
k] < 0.2*
radius && res_distances[
k] > 0.0 )
217 in->Set(TO_ERASE,
true);
232 tetgenio in,
out, in2, outnew;
237 in.numberofpoints = ThisModelPart.
Nodes().size();
238 in.pointlist =
new REAL[in.numberofpoints * 3];
241 ModelPart::NodesContainerType::iterator nodes_begin = ThisModelPart.
NodesBegin();
244 for(
unsigned int i = 0;
i<ThisModelPart.
Nodes().size();
i++)
246 (nodes_begin +
i)->SetId(
i+1);
251 for(
unsigned int i = 0;
i<ThisModelPart.
Nodes().size();
i++)
254 in.pointlist[base] = (nodes_begin +
i)->
X();
255 in.pointlist[base+1] = (nodes_begin +
i)->
Y();
256 in.pointlist[base+2] = (nodes_begin +
i)->
Z();
258 char tetgen_options[] =
"SQJ";
261 tetrahedralize(tetgen_options, &in, &
out);
264 double first_part_time = auxiliary.elapsed();
265 std::cout <<
"mesh generation time = " << first_part_time << std::endl;
268 int el_number =
out.numberoftetrahedra;
270 boost::timer alpha_shape_time;
272 std::vector<int> preserved_list(el_number);
275 int number_of_preserved_elems = 0;
276 for(
int el = 0;
el< el_number;
el++)
281 point_base = (
out.tetrahedronlist[base] - 1)*3;
282 x1[0] =
out.pointlist[point_base];
283 x1[1] =
out.pointlist[point_base+1];
284 x1[2] =
out.pointlist[point_base+2];
286 point_base = (
out.tetrahedronlist[base+1] - 1)*3;
287 x2[0] =
out.pointlist[point_base];
288 x2[1] =
out.pointlist[point_base+1];
289 x2[2] =
out.pointlist[point_base+2];
291 point_base = (
out.tetrahedronlist[base+2] - 1)*3;
292 x3[0] =
out.pointlist[point_base];
293 x3[1] =
out.pointlist[point_base+1];
294 x3[2] =
out.pointlist[point_base+2];
296 point_base = (
out.tetrahedronlist[base+3] - 1)*3;
297 x4[0] =
out.pointlist[point_base];
298 x4[1] =
out.pointlist[point_base+1];
299 x4[2] =
out.pointlist[point_base+2];
303 double geometrical_hmin, geometrical_hmax;
306 CalculateElementData(
x1,
x2,x3,x4,vol,
xc,
radius,geometrical_hmin,geometrical_hmax);
309 double prescribed_h = (nodes_begin +
out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(NODAL_H);
310 prescribed_h += (nodes_begin +
out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(NODAL_H);
311 prescribed_h += (nodes_begin +
out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(NODAL_H);
312 prescribed_h += (nodes_begin +
out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(NODAL_H);
313 prescribed_h *= 0.25;
316 int nb =
int( (nodes_begin +
out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
317 nb +=
int( (nodes_begin +
out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
318 nb +=
int( (nodes_begin +
out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
319 nb +=
int((nodes_begin +
out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_STRUCTURE) );
322 int nfs =
int( (nodes_begin +
out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
323 nfs +=
int( (nodes_begin +
out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
324 nfs +=
int( (nodes_begin +
out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
325 nfs +=
int((nodes_begin +
out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_FREE_SURFACE) );
328 int nfluid =
int( (nodes_begin +
out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_FLUID) );
329 nfluid +=
int( (nodes_begin +
out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_FLUID) );
330 nfluid +=
int( (nodes_begin +
out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_FLUID) );
331 nfluid +=
int((nodes_begin +
out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_FLUID) );
336 n_lag =
int( (nodes_begin +
out.tetrahedronlist[base]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
337 n_lag +=
int( (nodes_begin +
out.tetrahedronlist[base+1]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
338 n_lag +=
int( (nodes_begin +
out.tetrahedronlist[base+2]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
339 n_lag +=
int((nodes_begin +
out.tetrahedronlist[base+3]-1)->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET) );
366 preserved_list[
el] =
false;
372 if (n_lag >= 2 && nb==3 )
374 if(
radius < prescribed_h * 1.5*alpha_param )
377 preserved_list[
el] =
true;
378 number_of_preserved_elems += 1;
382 preserved_list[
el] =
false;
388 if(
radius < prescribed_h * alpha_param )
391 preserved_list[
el] =
true;
392 number_of_preserved_elems += 1;
396 preserved_list[
el] =
false;
402 else if (nfluid < 3.9)
404 if(
radius < prescribed_h * alpha_param )
407 preserved_list[
el] =
true;
408 number_of_preserved_elems += 1;
412 preserved_list[
el] =
false;
419 if(
radius < prescribed_h * alpha_param * 5.0 )
422 preserved_list[
el] =
true;
423 number_of_preserved_elems += 1;
428 preserved_list[
el] =
false;
437 std::cout <<
"time for passing alpha shape" << alpha_shape_time.elapsed() << std::endl;
447 in2.numberofpoints = ThisModelPart.
Nodes().size();
448 in2.pointlist =
new REAL[in2.numberofpoints * 3];
450 for(
unsigned int i = 0;
i<ThisModelPart.
Nodes().size();
i++)
453 in2.pointlist[base] = (nodes_begin +
i)->
X();
454 in2.pointlist[base+1] = (nodes_begin +
i)->
Y();
455 in2.pointlist[base+2] = (nodes_begin +
i)->
Z();
457 std::cout <<
"qui" << std::endl;
458 in2.numberoftetrahedra = number_of_preserved_elems;
459 in2.tetrahedronlist =
new int[in2.numberoftetrahedra * 4];
460 in2.tetrahedronvolumelist =
new double[in2.numberoftetrahedra];
463 for(
int el = 0;
el< el_number;
el++)
465 if( preserved_list[
el] )
470 in2.tetrahedronlist[new_base] =
out.tetrahedronlist[old_base];
471 in2.tetrahedronlist[new_base+1] =
out.tetrahedronlist[old_base+1];
472 in2.tetrahedronlist[new_base+2] =
out.tetrahedronlist[old_base+2];
473 in2.tetrahedronlist[new_base+3] =
out.tetrahedronlist[old_base+3];
476 double prescribed_h = (nodes_begin +
out.tetrahedronlist[old_base]-1)->FastGetSolutionStepValue(NODAL_H);
477 prescribed_h += (nodes_begin +
out.tetrahedronlist[old_base+1]-1)->FastGetSolutionStepValue(NODAL_H);
478 prescribed_h += (nodes_begin +
out.tetrahedronlist[old_base+2]-1)->FastGetSolutionStepValue(NODAL_H);
479 prescribed_h += (nodes_begin +
out.tetrahedronlist[old_base+3]-1)->FastGetSolutionStepValue(NODAL_H);
480 prescribed_h *= 0.25;
487 in2.tetrahedronvolumelist[
counter] = 2000.0*0.217*prescribed_h*prescribed_h*prescribed_h;
522 boost::timer mesh_recreation_time;
534 char mesh_regen_opts[] =
"rQJYq1.4/20nS";
535 tetrahedralize(mesh_regen_opts, &in2, &outnew);
540 char mesh_regen_opts[] =
"rQJYnS";
541 tetrahedralize(mesh_regen_opts, &in2, &outnew);
547 std::cout <<
"mesh recreation time" << mesh_recreation_time.elapsed() << std::endl;
577 int n_points_before_refinement = in2.numberofpoints;
579 if (outnew.numberofpoints>n_points_before_refinement)
581 for(
int i = n_points_before_refinement;
i<outnew.numberofpoints;
i++)
585 double&
x= outnew.pointlist[base];
586 double&
y= outnew.pointlist[base+1];
587 double&
z= outnew.pointlist[base+2];
593 list_of_new_nodes.push_back( pnode );
597 for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
602 (p_new_dof)->FreeDof();
608 std::cout <<
"During refinement we added " << outnew.numberofpoints-n_points_before_refinement<<
"nodes " <<std::endl;
650 if(outnew.numberofpoints-n_points_before_refinement > 0)
652 kd_tree nodes_tree2(list_of_new_nodes.begin(),list_of_new_nodes.end(),bucket_size);
656 for(
int el = 0;
el< in2.numberoftetrahedra;
el++)
663 point_base = (in2.tetrahedronlist[base] - 1)*3;
664 x1[0] = in2.pointlist[point_base];
665 x1[1] = in2.pointlist[point_base+1];
666 x1[2] = in2.pointlist[point_base+2];
668 point_base = (in2.tetrahedronlist[base+1] - 1)*3;
669 x2[0] = in2.pointlist[point_base];
670 x2[1] = in2.pointlist[point_base+1];
671 x2[2] = in2.pointlist[point_base+2];
673 point_base = (in2.tetrahedronlist[base+2] - 1)*3;
674 x3[0] = in2.pointlist[point_base];
675 x3[1] = in2.pointlist[point_base+1];
676 x3[2] = in2.pointlist[point_base+2];
678 point_base = (in2.tetrahedronlist[base+3] - 1)*3;
679 x4[0] = in2.pointlist[point_base];
680 x4[1] = in2.pointlist[point_base+1];
681 x4[2] = in2.pointlist[point_base+2];
685 double geometrical_hmin, geometrical_hmax;
694 CalculateElementData(
x1,
x2,x3,x4,vol,
xc,
radius,geometrical_hmin,geometrical_hmax);
698 std::size_t number_of_points_in_radius;
699 work_point.
X() =
xc[0]; work_point.
Y() =
xc[1]; work_point.
Z() =
xc[2];
701 number_of_points_in_radius = nodes_tree2.SearchInRadius(work_point,
radius, results.begin(),results_distances.begin(), max_results);
709 *( (nodes_begin + in2.tetrahedronlist[base]-1).base() ),
710 *( (nodes_begin + in2.tetrahedronlist[base+1]-1).base() ),
711 *( (nodes_begin + in2.tetrahedronlist[base+2]-1).base() ),
712 *( (nodes_begin + in2.tetrahedronlist[base+3]-1).base() )
717 for(
auto it=results.begin(); it!=results.begin() + number_of_points_in_radius; it++)
719 bool is_inside=
false;
721 is_inside = CalculatePosition(
x1[0],
x1[1],
x1[2],
725 (*it)->X(),(*it)->Y(), (*it)->Z(),
N);
727 if(is_inside ==
true)
730 Interpolate( geom,
N, step_data_size, *(it) );
742 for(
auto & list_of_new_node : list_of_new_nodes)
744 const array_1d<double,3>& disp = list_of_new_node->FastGetSolutionStepValue(DISPLACEMENT);
745 list_of_new_node->X0() = list_of_new_node->X() - disp[0];
746 list_of_new_node->Y0() = list_of_new_node->Y() - disp[1];
747 list_of_new_node->Z0() = list_of_new_node->Z() - disp[2];
755 boost::timer adding_elems;
759 (ThisModelPart.
Elements()).reserve(outnew.numberoftetrahedra);
761 for(
int iii = 0; iii< outnew.numberoftetrahedra; iii++)
766 *( (nodes_begin + outnew.tetrahedronlist[base]-1).base() ),
767 *( (nodes_begin + outnew.tetrahedronlist[base+1]-1).base() ),
768 *( (nodes_begin + outnew.tetrahedronlist[base+2]-1).base() ),
769 *( (nodes_begin + outnew.tetrahedronlist[base+3]-1).base() )
775 if( *(ModelNodes).find( outnew.tetrahedronlist[base]).base() == *(ThisModelPart.
Nodes().end()).base() )
777 if( *(ModelNodes).find( outnew.tetrahedronlist[base+1]).base() == *(ThisModelPart.
Nodes().end()).base() )
779 if( *(ModelNodes).find( outnew.tetrahedronlist[base+2]).base() == *(ThisModelPart.
Nodes().end()).base() )
781 if( *(ModelNodes).find( outnew.tetrahedronlist[base+3]).base() == *(ThisModelPart.
Nodes().end()).base() )
785 Element::Pointer p_element = rReferenceElement.
Create(
id, geom, properties);
786 (ThisModelPart.
Elements()).push_back(p_element);
789 std::cout <<
"time for adding elems" << adding_elems.elapsed() << std::endl;;
792 boost::timer adding_neighb;
794 ModelPart::ElementsContainerType::const_iterator el_begin = ThisModelPart.
ElementsBegin();
795 for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.
ElementsBegin();
799 int base = ( iii->Id() - 1 )*4;
801 (iii->GetValue(NEIGHBOUR_ELEMENTS)).resize(4);
804 for(
int i = 0;
i<4;
i++)
806 int index = outnew.neighborlist[base+
i];
810 neighb(
i) = Element::WeakPointer();
813 std::cout <<
"time for adding neigbours" << adding_neighb.elapsed() << std::endl;;
823 for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.
NodesBegin(); in!=ThisModelPart.
NodesEnd(); in++)
825 in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
831 boost::timer adding_faces;
833 (ThisModelPart.
Conditions()).reserve(outnew.numberoftrifaces );
836 for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.
ElementsBegin();
840 int base = ( iii->Id() - 1 )*4;
852 if( outnew.neighborlist[base] == -1)
854 CreateBoundaryFace(1, 2, 3, ThisModelPart, 0, *(iii.base()), rReferenceBoundaryCondition, properties );
857 if( outnew.neighborlist[base+1] == -1)
859 CreateBoundaryFace(0,3,2, ThisModelPart, 1, *(iii.base()), rReferenceBoundaryCondition, properties );
861 if( outnew.neighborlist[base+2] == -1)
864 CreateBoundaryFace(0,1,3, ThisModelPart, 2, *(iii.base()), rReferenceBoundaryCondition, properties );
866 if( outnew.neighborlist[base+3] == -1)
869 CreateBoundaryFace(0,2,1, ThisModelPart, 3, *(iii.base()), rReferenceBoundaryCondition, properties );
873 outnew.deinitialize();
875 std::cout <<
"time for adding faces" << adding_faces.elapsed() << std::endl;;
900 double second_part_time = auxiliary.elapsed();
901 std::cout <<
"second part time = " << second_part_time - first_part_time << std::endl;
922 std::string
Info()
const override{
return "";}
983 boost::numeric::ublas::bounded_matrix<double,3,3> mJ;
984 boost::numeric::ublas::bounded_matrix<double,3,3> mJinv;
989 void CreateBoundaryFace(
const int& i1,
const int& i2,
const int& i3,
ModelPart& ThisModelPart,
const int& outer_node_id, Element::Pointer origin_element,
Condition const& rReferenceBoundaryCondition, Properties::Pointer properties)
995 geom[i1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
996 geom[i2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
997 geom[i3].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
1002 temp.push_back(geom(i1));
1003 temp.push_back(geom(i2));
1004 temp.push_back(geom(i3));
1007 int id = (origin_element->Id()-1)*4;
1008 Condition::Pointer p_cond = rReferenceBoundaryCondition.
Create(
id,
temp, properties);
1011 (p_cond->GetValue(NEIGHBOUR_NODES)).clear();
1012 (p_cond->GetValue(NEIGHBOUR_NODES)).push_back( Node::WeakPointer( geom(outer_node_id) ) );
1013 (p_cond->GetValue(NEIGHBOUR_ELEMENTS)).clear();
1014 (p_cond->GetValue(NEIGHBOUR_ELEMENTS)).push_back( Element::WeakPointer( origin_element ) );
1015 ThisModelPart.
Conditions().push_back(p_cond);
1022 void Interpolate( Tetrahedra3D4<Node >& geom,
const array_1d<double,4>&
N,
1023 unsigned int step_data_size,
1024 Node::Pointer pnode)
1026 unsigned int buffer_size = pnode->GetBufferSize();
1033 double* step_data = (pnode)->SolutionStepData().Data(
step);
1036 double* node0_data = geom[0].SolutionStepData().Data(
step);
1037 double* node1_data = geom[1].SolutionStepData().Data(
step);
1038 double* node2_data = geom[2].SolutionStepData().Data(
step);
1039 double* node3_data = geom[3].SolutionStepData().Data(
step);
1042 for(
unsigned int j= 0;
j<step_data_size;
j++)
1045 step_data[
j] =
N[0]*node0_data[
j] +
N[1]*node1_data[
j] +
N[2]*node2_data[
j] +
N[3]*node3_data[
j];
1055 pnode->FastGetSolutionStepValue(IS_BOUNDARY)=0.0;
1056 pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1057 pnode->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET)=0.0;
1058 pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=0.0;
1059 pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1062 inline double CalculateVol(
const double x0,
const double y0,
const double z0,
1063 const double x1,
const double y1,
const double z1,
1064 const double x2,
const double y2,
const double z2,
1065 const double x3,
const double y3,
const double z3
1068 double x10 =
x1 - x0;
1069 double y10 = y1 - y0;
1070 double z10 = z1 - z0;
1072 double x20 =
x2 - x0;
1073 double y20 = y2 - y0;
1074 double z20 = z2 - z0;
1076 double x30 = x3 - x0;
1077 double y30 = y3 - y0;
1078 double z30 = z3 - z0;
1080 double detJ = x10 * y20 * z30 - x10 * y30 * z20 + y10 * z20 * x30 - y10 * x20 * z30 + z10 * x20 * y30 - z10 * y20 * x30;
1081 return detJ*0.1666666666666666666667;
1086 inline bool CalculatePosition(
const double x0,
const double y0,
const double z0,
1087 const double x1,
const double y1,
const double z1,
1088 const double x2,
const double y2,
const double z2,
1089 const double x3,
const double y3,
const double z3,
1090 const double xc,
const double yc,
const double zc,
1091 array_1d<double,4>&
N
1096 double vol = CalculateVol(x0,y0,z0,
x1,y1,z1,
x2,y2,z2,x3,y3,z3);
1098 double inv_vol = 0.0;
1102 KRATOS_WATCH(
"WARNING!!!!!!!!!! YOU HAVE A SLIVER HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
1107 inv_vol = 1.0 / vol;
1108 N[0] = CalculateVol(
x1,y1,z1,x3,y3,z3,
x2,y2,z2,
xc,
yc,zc) * inv_vol;
1109 N[1] = CalculateVol(x3,y3,z3,x0,y0,z0,
x2,y2,z2,
xc,
yc,zc) * inv_vol;
1110 N[2] = CalculateVol(x3,y3,z3,
x1,y1,z1,x0,y0,z0,
xc,
yc,zc) * inv_vol;
1111 N[3] = CalculateVol(x0,y0,z0,
x1,y1,z1,
x2,y2,z2,
xc,
yc,zc) * inv_vol;
1114 if(
N[0] >= 0.0 &&
N[1] >= 0.0 &&
N[2] >= 0.0 &&
N[3] >=0.0 &&
N[0] <= 1.0 &&
N[1] <= 1.0 &&
N[2] <= 1.0 &&
N[3] <=1.0)
1127 inline void CalculateCenterAndSearchRadius(
const double x0,
const double y0,
const double z0,
1128 const double x1,
const double y1,
const double z1,
1129 const double x2,
const double y2,
const double z2,
1130 const double x3,
const double y3,
const double z3,
1131 double&
xc,
double&
yc,
double& zc,
double&
R
1135 yc = 0.25*(y0+y1+y2+y3);
1136 zc = 0.25*(z0+z1+z2+z3);
1138 double R1 = (
xc-x0)*(
xc-x0) + (
yc-y0)*(
yc-y0) + (zc-z0)*(zc-z0);
1139 double R2 = (
xc-
x1)*(
xc-
x1) + (
yc-y1)*(
yc-y1) + (zc-z1)*(zc-z1);
1140 double R3 = (
xc-
x2)*(
xc-
x2) + (
yc-y2)*(
yc-y2) + (zc-z2)*(zc-z2);
1141 double R4 = (
xc-x3)*(
xc-x3) + (
yc-y3)*(
yc-y3) + (zc-z3)*(zc-z3);
1152 void CalculateElementData(
const array_1d<double,3>& c1,
1153 const array_1d<double,3>& c2,
1154 const array_1d<double,3>& c3,
1155 const array_1d<double,3>& c4,
1157 array_1d<double,3>&
xc,
1164 const double x0 = c1[0];
const double y0 = c1[1];
const double z0 = c1[2];
1165 const double x1 = c2[0];
const double y1 = c2[1];
const double z1 = c2[2];
1166 const double x2 = c3[0];
const double y2 = c3[1];
const double z2 = c3[2];
1167 const double x3 = c4[0];
const double y3 = c4[1];
const double z3 = c4[2];
1183 mJ(0,0) =
x1-x0; mJ(0,1) = y1-y0; mJ(0,2) = z1-z0;
1184 mJ(1,0) =
x2-x0; mJ(1,1) = y2-y0; mJ(1,2) = z2-z0;
1185 mJ(2,0) = x3-x0; mJ(2,1) = y3-y0; mJ(2,2) = z3-z0;
1190 mJinv(0,0) = mJ(1,1)*mJ(2,2) - mJ(1,2)*mJ(2,1);
1191 mJinv(1,0) = -mJ(1,0)*mJ(2,2) + mJ(1,2)*mJ(2,0);
1192 mJinv(2,0) = mJ(1,0)*mJ(2,1) - mJ(1,1)*mJ(2,0);
1194 mJinv(0,1) = -mJ(0,1)*mJ(2,2) + mJ(0,2)*mJ(2,1);
1195 mJinv(1,1) = mJ(0,0)*mJ(2,2) - mJ(0,2)*mJ(2,0);
1196 mJinv(2,1) = -mJ(0,0)*mJ(2,1) + mJ(0,1)*mJ(2,0);
1198 mJinv(0,2) = mJ(0,1)*mJ(1,2) - mJ(0,2)*mJ(1,1);
1199 mJinv(1,2) = -mJ(0,0)*mJ(1,2) + mJ(0,2)*mJ(1,0);
1200 mJinv(2,2) = mJ(0,0)*mJ(1,1) - mJ(0,1)*mJ(1,0);
1204 double detJ = mJ(0,0)*mJinv(0,0)
1205 + mJ(0,1)*mJinv(1,0)
1206 + mJ(0,2)*mJinv(2,0);
1208 volume = detJ * 0.16666666667;
1211 if(
volume < 1
e-3 * hmax*hmax*hmax)
1215 xc[0]=0.0;
xc[1]=0.0;
xc[2]=0.0;
1223 double x0_2 = x0*x0 + y0*y0 + z0*z0;
1224 double x1_2 =
x1*
x1 + y1*y1 + z1*z1;
1225 double x2_2 =
x2*
x2 + y2*y2 + z2*z2;
1226 double x3_2 = x3*x3 + y3*y3 + z3*z3;
1233 mRhs[0] = 0.5*(x1_2 - x0_2);
1234 mRhs[1] = 0.5*(x2_2 - x0_2);
1235 mRhs[2] = 0.5*(x3_2 - x0_2);
1250 template<
class T, std::
size_t dim >
1251 class PointDistance{
1253 double operator()( T
const& p1, T
const& p2 ){
1255 for( std::size_t
i = 0 ;
i <
dim ;
i++){
1256 double tmp = p1[
i] - p2[
i];
1263 template<
class T, std::
size_t dim >
1264 class DistanceCalculator{
1266 double operator()( T
const& p1, T
const& p2 ){
1268 for( std::size_t
i = 0 ;
i <
dim ;
i++){
1269 double tmp = p1[
i] - p2[
i];
1327 rOStream << std::endl;
Short class definition.
Definition: bucket.h:57
Base class for all Conditions.
Definition: condition.h:59
virtual Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const
It creates a new condition pointer.
Definition: condition.h:205
Dof represents a degree of freedom (DoF).
Definition: dof.h:86
Base class for all Elements.
Definition: element.h:60
virtual Pointer Create(IndexType NewId, NodesArrayType const &ThisNodes, PropertiesType::Pointer pProperties) const
It creates a new element pointer.
Definition: element.h:202
This process removes the entities from a model part with the flag TO_ERASE.
Definition: entity_erase_process.h:70
void Execute() override
Execute method is used to execute the Process algorithms.
Geometry base class.
Definition: geometry.h:71
This class is a wrapper for a pointer to a data that is located in a different rank.
Definition: global_pointer.h:44
Definition: amatrix_interface.h:497
Definition: amatrix_interface.h:530
PropertiesType::Pointer pGetProperties(IndexType PropertiesId)
Definition: mesh.h:394
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
NodeType::Pointer CreateNewNode(int Id, double x, double y, double z, VariablesList::Pointer pNewVariablesList, IndexType ThisIndex=0)
Definition: model_part.cpp:270
SizeType GetNodalSolutionStepDataSize()
Definition: model_part.h:571
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
ConditionsContainerType & Conditions(IndexType ThisIndex=0)
Definition: model_part.h:1381
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
ElementIterator ElementsEnd(IndexType ThisIndex=0)
Definition: model_part.h:1179
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
MeshType & GetMesh(IndexType ThisIndex=0)
Definition: model_part.h:1791
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
VariablesList & GetNodalSolutionStepVariablesList()
Definition: model_part.h:549
This class defines the node.
Definition: node.h:65
std::vector< std::unique_ptr< Dof< double > >> DofsContainerType
The DoF container type definition.
Definition: node.h:92
double Y() const
Definition: point.h:187
double Z() const
Definition: point.h:193
double X() const
Definition: point.h:181
PointerVector is a container like stl vector but using a vector to store pointers to its data.
Definition: pointer_vector.h:72
Short class definition.
Definition: tetgen_pfem_refine.h:67
Short class definition.
Definition: tetgen_pfem_refine_vms.h:68
void PrintInfo(std::ostream &rOStream) const override
Print information about this object.
Definition: tetgen_pfem_refine_vms.h:925
~TetGenPfemModelerVms() override=default
Destructor.
TetGenPfemModelerVms()
Default constructor.
Definition: tetgen_pfem_refine_vms.h:81
void ReGenerateMesh(ModelPart &ThisModelPart, Element const &rReferenceElement, Condition const &rReferenceBoundaryCondition, EntitiesEraseProcess< Node > &node_erase, bool rem_nodes=true, bool add_nodes=true, double alpha_param=1.4, double h_factor=0.5)
Definition: tetgen_pfem_refine_vms.h:103
void PrintData(std::ostream &rOStream) const override
Print object's data.
Definition: tetgen_pfem_refine_vms.h:928
KRATOS_CLASS_POINTER_DEFINITION(TetGenPfemModelerVms)
Pointer definition of TetGenPfemModeler.
std::string Info() const override
Turn back information as a string.
Definition: tetgen_pfem_refine_vms.h:922
A four node tetrahedra geometry with linear shape functions.
Definition: tetrahedra_3d_4.h:59
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
A three node 3D triangle geometry with linear shape functions.
Definition: triangle_3d_3.h:77
bool Has(const VariableData &rThisVariable) const
Definition: variables_list.h:372
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_WATCH(variable)
Definition: define.h:806
#define KRATOS_TRY
Definition: define.h:109
#define REAL
Definition: delaunator_utilities.cpp:16
z
Definition: GenerateWind.py:163
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
DistanceVector::iterator DistanceIterator
Definition: internal_variables_interpolation_process.h:56
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
PointVector::iterator PointIterator
Definition: internal_variables_interpolation_process.h:54
std::ostream & operator<<(std::ostream &rOStream, const LinearMasterSlaveConstraint &rThis)
output stream function
Definition: linear_master_slave_constraint.h:432
std::vector< double > DistanceVector
Definition: internal_variables_interpolation_process.h:55
Point PointType
Geometric definitions.
Definition: mortar_classes.h:36
float dist
Definition: edgebased_PureConvection.py:89
int step
Definition: face_heat.py:88
y
Other simbols definition.
Definition: generate_axisymmetric_navier_stokes_element.py:54
res
Definition: generate_convection_diffusion_explicit_element.py:211
x2
Definition: generate_frictional_mortar_condition.py:122
x1
Definition: generate_frictional_mortar_condition.py:121
tuple tmp
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:98
out
Definition: isotropic_damage_automatic_differentiation.py:200
R
Definition: isotropic_damage_automatic_differentiation.py:172
float radius
Definition: mesh_to_mdpa_converter.py:18
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
el
Definition: read_stl.py:25
float temp
Definition: rotating_cone.py:85
float xc
Definition: rotating_cone.py:77
float yc
Definition: rotating_cone.py:78
int counter
Definition: script_THERMAL_CORRECT.py:218
add_nodes
Definition: script.py:96
N
Definition: sensitivityMatrix.py:29
x
Definition: sensitivityMatrix.py:49
int dim
Definition: sensitivityMatrix.py:25
volume
Definition: sp_statistics.py:15
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31