13 #if !defined(KRATOS_TRIGEN_GLASS_MODELER_H_INCLUDED )
14 #define KRATOS_TRIGEN_GLASS_MODELER_H_INCLUDED
33 void triangulate(
char *,
struct triangulateio *,
struct triangulateio *,
struct triangulateio *);
107 Element const& rReferenceElement,
108 Condition const& rReferenceBoundaryCondition,
110 double my_alpha = 1.4,
double h_factor=0.5)
114 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_FREE_SURFACE)==
false )
115 KRATOS_THROW_ERROR(std::logic_error,
"Add ----IS_FREE_SURFACE---- variable!!!!!! ERROR",
"");
116 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_STRUCTURE)==
false )
117 KRATOS_THROW_ERROR(std::logic_error,
"Add ----IS_STRUCTURE---- variable!!!!!! ERROR",
"");
118 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_BOUNDARY)==
false )
119 KRATOS_THROW_ERROR(std::logic_error,
"Add ----IS_BOUNDARY---- variable!!!!!! ERROR",
"");
120 if (ThisModelPart.
NodesBegin()->SolutionStepsDataHas(IS_FLUID)==
false )
124 const auto inital_time = std::chrono::steady_clock::now();
133 int id = (ThisModelPart.
Nodes().end() - 1)->Id() + 1;
137 if (ic->GetGeometry().size()==2)
139 unsigned int n_fs=ic->GetGeometry()[0].FastGetSolutionStepValue(IS_FREE_SURFACE);
140 n_fs+=ic->GetGeometry()[1].FastGetSolutionStepValue(IS_FREE_SURFACE);
143 if (n_fs==ic->GetGeometry().size())
145 double x0=ic->GetGeometry()[0].X();
double y0=ic->GetGeometry()[0].Y();
146 double x1=ic->GetGeometry()[1].X();
double y1=ic->GetGeometry()[1].Y();
149 double edge=sqrt((x0-
x1)*(x0-
x1)+(y0-y1)*(y0-y1));
151 double nodal_h=ic->GetGeometry()[0].FastGetSolutionStepValue(NODAL_H);
152 nodal_h+=ic->GetGeometry()[1].FastGetSolutionStepValue(NODAL_H);
165 double x=0.5*(x0+
x1);
double y=0.5*(y0+y1);
double z=0.0;
174 for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
179 (p_new_dof)->FreeDof();
182 InterpolateOnEdge(geom, step_data_size, pnode);
184 pnode->X0() = pnode->X() - disp[0];
185 pnode->Y0() = pnode->Y() - disp[1];
210 unsigned int bucket_size = 20;
213 unsigned int max_results = 100;
218 Node work_point(0,0.0,0.0,0.0);
224 list_of_nodes.reserve(ThisModelPart.
Nodes().size());
225 for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.
NodesBegin() ; i_node != ThisModelPart.
NodesEnd() ; i_node++)
227 (list_of_nodes).push_back(*(i_node.base()));
230 KdtreeType nodes_tree1(list_of_nodes.begin(),list_of_nodes.end(), bucket_size);
232 RemoveCloseNodes(ThisModelPart, nodes_tree1, node_erase, h_factor);
240 struct triangulateio in_mid;
241 struct triangulateio out_mid;
242 struct triangulateio vorout_mid;
244 initialize_triangulateio(in_mid);
245 initialize_triangulateio(out_mid);
246 initialize_triangulateio(vorout_mid);
250 in_mid.numberofpoints = ThisModelPart.
Nodes().size();
251 in_mid.pointlist = (
REAL *) malloc(in_mid.numberofpoints * 2 *
sizeof(
REAL));
254 ModelPart::NodesContainerType::iterator nodes_begin = ThisModelPart.
NodesBegin();
255 for(
unsigned int i = 0;
i<ThisModelPart.
Nodes().size();
i++)
261 (nodes_begin +
i)->SetId(
i+1);
264 in_mid.pointlist[base] = (nodes_begin +
i)->
X();
265 in_mid.pointlist[base+1] = (nodes_begin +
i)->
Y();
267 auto& node_dofs = (nodes_begin +
i)->GetDofs();
268 for(
auto iii = node_dofs.begin(); iii != node_dofs.end(); iii++)
270 (**iii).SetEquationId(
i+1);
281 char options1[] =
"Pne";
282 triangulate(options1, &in_mid, &out_mid, &vorout_mid);
286 unsigned int el_number=out_mid.numberoftriangles;
289 std::vector<int> preserved_list1(el_number);
290 preserved_list1.resize(el_number);
292 int number_of_preserved_elems= PassAlphaShape(ThisModelPart, preserved_list1, el_number, my_alpha, out_mid);
296 clean_triangulateio(in_mid);
297 clean_triangulateio(vorout_mid);
301 struct triangulateio in2;
302 struct triangulateio out2;
303 struct triangulateio vorout2;
305 initialize_triangulateio(in2);
306 initialize_triangulateio(out2);
307 initialize_triangulateio(vorout2);
310 in2.numberofpoints = ThisModelPart.
Nodes().size();
311 in2.pointlist = (
REAL *) malloc(in2.numberofpoints * 2 *
sizeof(
REAL));
314 for(
unsigned int i = 0;
i<ThisModelPart.
Nodes().size();
i++)
317 in2.pointlist[base] = (nodes_begin +
i)->
X();
318 in2.pointlist[base+1] = (nodes_begin +
i)->
Y();
320 in2.numberoftriangles=number_of_preserved_elems;
322 in2.trianglelist = (
int *) malloc(in2.numberoftriangles * 3 *
sizeof(
int));
323 in2.trianglearealist = (
REAL *) malloc(in2.numberoftriangles *
sizeof(
REAL));
330 for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.
NodesBegin() ; i_node != ThisModelPart.
NodesEnd() ; i_node++)
332 if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
335 double& val=i_node->FastGetSolutionStepValue(NODAL_H);
343 for (
unsigned int el=0;
el<el_number;
el++)
345 if(
static_cast<bool>(preserved_list1[
el]) == true )
351 in2.trianglelist[new_base] = out_mid.trianglelist[old_base];
352 in2.trianglelist[new_base+1] = out_mid.trianglelist[old_base+1];
353 in2.trianglelist[new_base+2] = out_mid.trianglelist[old_base+2];
356 double prescribed_h = (nodes_begin + out_mid.trianglelist[old_base]-1)->FastGetSolutionStepValue(NODAL_H);
357 prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+1]-1)->FastGetSolutionStepValue(NODAL_H);
358 prescribed_h += (nodes_begin + out_mid.trianglelist[old_base+2]-1)->FastGetSolutionStepValue(NODAL_H);
359 prescribed_h *= 0.3333;
361 in2.trianglearealist[
counter] = 0.5*(1.5*prescribed_h*1.5*prescribed_h);
367 for(ModelPart::NodesContainerType::iterator i_node = ThisModelPart.
NodesBegin() ; i_node != ThisModelPart.
NodesEnd() ; i_node++)
369 if (i_node->FastGetSolutionStepValue(IS_FREE_SURFACE)!=0)
371 double& nodal_h=i_node->FastGetSolutionStepValue(NODAL_H);
377 clean_triangulateio(out_mid);
386 char mesh_regen_opts[] =
"YJq1.4arn";
387 triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
392 char mesh_regen_opts[] =
"YJrn";
393 triangulate(mesh_regen_opts, &in2, &out2, &vorout2);
401 typedef std::vector<PointType::Pointer>
PointVector;
421 int n_points_before_refinement = in2.numberofpoints;
423 if (out2.numberofpoints > n_points_before_refinement )
425 for(
int i = n_points_before_refinement;
i<out2.numberofpoints;
i++)
429 double&
x= out2.pointlist[base];
430 double&
y= out2.pointlist[base+1];
434 pnode->SetBufferSize(ThisModelPart.
NodesBegin()->GetBufferSize() );
436 list_of_new_nodes.push_back( pnode );
440 for(Node::DofsContainerType::iterator iii = reference_dofs.begin(); iii != reference_dofs.end(); iii++)
445 (p_new_dof)->FreeDof();
453 std::cout <<
"During refinement we added " << out2.numberofpoints-n_points_before_refinement<<
" nodes " <<std::endl;
468 unsigned int MaximumNumberOfResults = list_of_new_nodes.size();
474 if(out2.numberofpoints-n_points_before_refinement > 0)
477 kd_tree nodes_tree2(list_of_new_nodes.begin(),list_of_new_nodes.end(),bucket_size);
480 for(
int el = 0;
el< in2.numberoftriangles;
el++)
484 point_base = (in2.trianglelist[base] - 1)*2;
485 x1[0] = in2.pointlist[point_base];
486 x1[1] = in2.pointlist[point_base+1];
488 point_base = (in2.trianglelist[base+1] - 1)*2;
489 x2[0] = in2.pointlist[point_base];
490 x2[1] = in2.pointlist[point_base+1];
492 point_base = (in2.trianglelist[base+2] - 1)*2;
493 x3[0] = in2.pointlist[point_base];
494 x3[1] = in2.pointlist[point_base+1];
499 CalculateCenterAndSearchRadius(
x1[0],
x1[1],
505 int number_of_points_in_radius;
508 work_point.
Z() = 0.0;
510 number_of_points_in_radius = nodes_tree2.SearchInRadius(work_point,
radius*1.01, Results.begin(),
511 ResultsDistances.begin(), MaximumNumberOfResults);
514 *( (nodes_begin + in2.trianglelist[base]-1).base() ),
515 *( (nodes_begin + in2.trianglelist[base+1]-1).base() ),
516 *( (nodes_begin + in2.trianglelist[base+2]-1).base() )
520 for(
PointIterator it_found = Results.begin(); it_found != Results.begin() + number_of_points_in_radius; it_found++)
522 bool is_inside =
false;
526 (*it_found)->X(),(*it_found)->Y(),
N);
529 if(is_inside ==
true)
541 for( PointVector::iterator it = list_of_new_nodes.begin(); it!=list_of_new_nodes.end(); it++)
544 (*it)->X0() = (*it)->X() - disp[0];
545 (*it)->Y0() = (*it)->Y() - disp[1];
556 (ThisModelPart.
Elements()).reserve(out2.numberoftriangles);
558 for(
int iii = 0; iii< out2.numberoftriangles; iii++)
563 *( (nodes_begin + out2.trianglelist[base]-1).base() ),
564 *( (nodes_begin + out2.trianglelist[base+1]-1).base() ),
565 *( (nodes_begin + out2.trianglelist[base+2]-1).base() )
571 if( *(ModelNodes).find( out2.trianglelist[base]).base() == *(ThisModelPart.
Nodes().end()).base() )
573 if( *(ModelNodes).find( out2.trianglelist[base+1]).base() == *(ThisModelPart.
Nodes().end()).base() )
575 if( *(ModelNodes).find( out2.trianglelist[base+2]).base() == *(ThisModelPart.
Nodes().end()).base() )
580 Element::Pointer p_element = rReferenceElement.
Create(
id, geom, properties);
581 (ThisModelPart.
Elements()).push_back(p_element);
587 ModelPart::ElementsContainerType::const_iterator el_begin = ThisModelPart.
ElementsBegin();
588 for(ModelPart::ElementsContainerType::const_iterator iii = ThisModelPart.
ElementsBegin();
592 int base = ( iii->Id() - 1 )*3;
594 (iii->GetValue(NEIGHBOUR_ELEMENTS)).resize(3);
596 for(
int i = 0;
i<3;
i++)
598 int index = out2.neighborlist[base+
i];
602 neighb(
i) = Element::WeakPointer();
606 IdentifyBoundary(ThisModelPart, rReferenceBoundaryCondition, properties, out2);
610 clean_triangulateio(in2);
612 clean_triangulateio(out2);
614 clean_triangulateio(vorout2);
615 KRATOS_WATCH(
"Finished remeshing with Trigen_GLASS_Refine")
636 virtual std::string
Info()
const
700 boost::numeric::ublas::bounded_matrix<double,2,2> mJ;
701 boost::numeric::ublas::bounded_matrix<double,2,2> mJinv;
715 unsigned int max_results = 100;
720 Node work_point(0,0.0,0.0,0.0);
722 unsigned int n_points_in_radius;
726 for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.
NodesBegin(); in != ThisModelPart.
NodesEnd(); in++)
728 radius=h_factor*in->FastGetSolutionStepValue(NODAL_H);
730 work_point[0]=in->X();
731 work_point[1]=in->Y();
732 work_point[2]=in->Z();
734 n_points_in_radius = nodes_tree1.SearchInRadius(work_point,
radius,
res.begin(),res_distances.begin(), max_results);
735 if (n_points_in_radius>1)
737 if (in->FastGetSolutionStepValue(IS_BOUNDARY)==0.0 && in->FastGetSolutionStepValue(IS_STRUCTURE)==0.0)
740 double erased_nodes = 0;
742 erased_nodes += in->Is(TO_ERASE);
744 if( erased_nodes < 1)
745 in->Set(TO_ERASE,
true);
749 else if ( (in)->FastGetSolutionStepValue(IS_STRUCTURE)!=1.0 && (in)->FastGetSolutionStepValue(FLAG_VARIABLE)!=1.0)
757 if ( (*i)->FastGetSolutionStepValue(IS_BOUNDARY,1)==1.0 && res_distances[
k] < 0.2*
radius && res_distances[
k] > 0.0 )
765 in->Set(TO_ERASE,
true);
779 int PassAlphaShape(
ModelPart& ThisModelPart,
std::vector<
int>& preserved_list1,
unsigned int & el_number,
double& my_alpha, struct triangulateio& out_mid)
787 array_1d<double,3>
x1,
x2,x3,
xc;
790 int number_of_preserved_elems=0;
793 for(
unsigned int el = 0;
el< el_number;
el++)
798 point_base = (out_mid.trianglelist[base] - 1)*2;
799 x1[0] = out_mid.pointlist[point_base];
800 x1[1] = out_mid.pointlist[point_base+1];
802 point_base = (out_mid.trianglelist[base+1] - 1)*2;
803 x2[0] = out_mid.pointlist[point_base];
804 x2[1] = out_mid.pointlist[point_base+1];
806 point_base = (out_mid.trianglelist[base+2] - 1)*2;
807 x3[0] = out_mid.pointlist[point_base];
808 x3[1] = out_mid.pointlist[point_base+1];
811 Geometry<Node >
temp;
813 temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base]).base() ) );
814 temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+1]).base() ) );
815 temp.push_back( *((ThisModelPart.Nodes()).find( out_mid.trianglelist[base+2]).base() ) );
817 int number_of_structure_nodes =
int(
temp[0].FastGetSolutionStepValue(IS_STRUCTURE) );
818 number_of_structure_nodes +=
int(
temp[1].FastGetSolutionStepValue(IS_STRUCTURE) );
819 number_of_structure_nodes +=
int(
temp[2].FastGetSolutionStepValue(IS_STRUCTURE) );
822 int nfs =
int(
temp[0].FastGetSolutionStepValue(IS_FREE_SURFACE) );
823 nfs +=
int(
temp[1].FastGetSolutionStepValue(IS_FREE_SURFACE) );
824 nfs +=
int(
temp[2].FastGetSolutionStepValue(IS_FREE_SURFACE) );
827 int nfluid =
int(
temp[0].FastGetSolutionStepValue(IS_FLUID) );
828 nfluid +=
int(
temp[1].FastGetSolutionStepValue(IS_FLUID) );
829 nfluid +=
int(
temp[2].FastGetSolutionStepValue(IS_FLUID) );
846 if(number_of_structure_nodes!=3)
850 if (nfs != 0 || nfluid != 3)
853 if( AlphaShape(my_alpha,
temp))
861 preserved_list1[
el] =
true;
862 number_of_preserved_elems += 1;
884 double bigger_alpha = my_alpha*10.0;
885 if( AlphaShape(bigger_alpha,
temp) && number_of_structure_nodes!=3)
888 preserved_list1[
el] =
true;
889 number_of_preserved_elems += 1;
896 preserved_list1[
el] =
false;
899 return number_of_preserved_elems;
906 void IdentifyBoundary(
ModelPart& ThisModelPart, Condition
const& rReferenceBoundaryCondition, Properties::Pointer& properties,
struct triangulateio& out2)
910 for(ModelPart::NodesContainerType::const_iterator in = ThisModelPart.NodesBegin(); in!=ThisModelPart.NodesEnd(); in++)
912 in->FastGetSolutionStepValue(IS_BOUNDARY) = 0;
915 ModelPart::ElementsContainerType::iterator elements_end = ThisModelPart.Elements().end();
917 ThisModelPart.Elements().Unique();
920 for(ModelPart::ElementsContainerType::iterator iii = ThisModelPart.ElementsBegin(); iii != ThisModelPart.ElementsEnd(); iii++)
922 int base = ( iii->Id() - 1 )*3;
924 ModelPart::ElementsContainerType::iterator el_neighb;
935 el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base] );
936 if( el_neighb == elements_end )
940 iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
941 iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
946 temp1.push_back(iii->GetGeometry()(1));
947 temp1.push_back(iii->GetGeometry()(2));
949 Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(
new Geometry< Node >(temp1) );
950 int id = (iii->Id()-1)*3;
951 Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(
id, temp1, properties);
952 ThisModelPart.Conditions().push_back(p_cond);
958 el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+1] );
960 if( el_neighb == elements_end )
963 iii->GetGeometry()[2].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
964 iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
969 temp1.push_back(iii->GetGeometry()(2));
970 temp1.push_back(iii->GetGeometry()(0));
972 Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(
new Geometry< Node >(temp1) );
973 int id = (iii->Id()-1)*3+1;
974 Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(
id, temp1, properties);
975 ThisModelPart.Conditions().push_back(p_cond);
982 el_neighb = (ThisModelPart.Elements()).find( out2.neighborlist[base+2] );
983 if( el_neighb == elements_end )
986 iii->GetGeometry()[0].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
987 iii->GetGeometry()[1].FastGetSolutionStepValue(IS_BOUNDARY) = 1;
992 temp1.push_back(iii->GetGeometry()(0));
993 temp1.push_back(iii->GetGeometry()(1));
995 Geometry< Node >::Pointer cond = Geometry< Node >::Pointer(
new Geometry< Node >(temp1) );
996 int id = (iii->Id()-1)*3+2;
997 Condition::Pointer p_cond = rReferenceBoundaryCondition.Create(
id, temp1, properties);
998 ThisModelPart.Conditions().push_back(p_cond);
1008 bool AlphaShape(
double alpha_param, Geometry<Node >& pgeom)
1014 double x0 = pgeom[0].X();
1015 double x1 = pgeom[1].X();
1016 double x2 = pgeom[2].X();
1018 double y0 = pgeom[0].Y();
1019 double y1 = pgeom[1].Y();
1020 double y2 = pgeom[2].Y();
1022 mJ(0,0)=2.0*(
x1-x0);
1023 mJ(0,1)=2.0*(y1-y0);
1024 mJ(1,0)=2.0*(
x2-x0);
1025 mJ(1,1)=2.0*(y2-y0);
1028 double detJ = mJ(0,0)*mJ(1,1)-mJ(0,1)*mJ(1,0);
1030 mJinv(0,0) = mJ(1,1);
1031 mJinv(0,1) = -mJ(0,1);
1032 mJinv(1,0) = -mJ(1,0);
1033 mJinv(1,1) = mJ(0,0);
1035 bounded_matrix<double,2,2> check;
1039 h = pgeom[0].FastGetSolutionStepValue(NODAL_H);
1040 h += pgeom[1].FastGetSolutionStepValue(NODAL_H);
1041 h += pgeom[2].FastGetSolutionStepValue(NODAL_H);
1049 pgeom[0].GetSolutionStepValue(IS_BOUNDARY) = 1;
1050 pgeom[1].GetSolutionStepValue(IS_BOUNDARY) = 1;
1051 pgeom[2].GetSolutionStepValue(IS_BOUNDARY) = 1;
1058 double x0_2 = x0*x0 + y0*y0;
1059 double x1_2 =
x1*
x1 + y1*y1;
1060 double x2_2 =
x2*
x2 + y2*y2;
1066 mRhs[0] = (x1_2 - x0_2);
1067 mRhs[1] = (x2_2 - x0_2);
1072 double radius = sqrt(pow(mC[0]-x0,2)+pow(mC[1]-y0,2));
1089 inline void CalculateCenterAndSearchRadius(
const double x0,
const double y0,
1090 const double x1,
const double y1,
1091 const double x2,
const double y2,
1092 double&
xc,
double&
yc,
double&
R
1095 xc = 0.3333333333333333333*(x0+
x1+
x2);
1096 yc = 0.3333333333333333333*(y0+y1+y2);
1098 double R1 = (
xc-x0)*(
xc-x0) + (
yc-y0)*(
yc-y0);
1109 inline double CalculateVol(
const double x0,
const double y0,
1110 const double x1,
const double y1,
1111 const double x2,
const double y2
1114 return 0.5*( (
x1-x0)*(y2-y0)- (y1-y0)*(
x2-x0) );
1118 const double x1,
const double y1,
1119 const double x2,
const double y2,
1120 const double xc,
const double yc,
1121 array_1d<double,3>&
N
1126 if(area < 0.000000000001)
1141 double upper_limit = 1.0+
tol;
1142 double lower_limit = -
tol;
1144 if(
N[0] >= lower_limit &&
N[1] >= lower_limit &&
N[2] >= lower_limit &&
N[0] <= upper_limit &&
N[1] <= upper_limit &&
N[2] <= upper_limit)
1151 void Interpolate( Triangle2D3<Node >& geom,
const array_1d<double,3>&
N,
1152 unsigned int step_data_size,
1153 Node::Pointer pnode)
1155 unsigned int buffer_size = pnode->GetBufferSize();
1163 double* step_data = (pnode)->SolutionStepData().Data(
step);
1166 double* node0_data = geom[0].SolutionStepData().Data(
step);
1167 double* node1_data = geom[1].SolutionStepData().Data(
step);
1168 double* node2_data = geom[2].SolutionStepData().Data(
step);
1171 for(
unsigned int j= 0;
j<step_data_size;
j++)
1174 step_data[
j] =
N[0]*node0_data[
j] +
N[1]*node1_data[
j] +
N[2]*node2_data[
j];
1179 if (
N[0]==0.0 &&
N[1]==0.0 &&
N[2]==0.0)
1180 KRATOS_THROW_ERROR(std::logic_error,
"SOMETHING's wrong with the added nodes!!!!!! ERROR",
"");
1191 pnode->FastGetSolutionStepValue(IS_BOUNDARY)=0.0;
1192 pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1193 pnode->FastGetSolutionStepValue(IS_INTERFACE)=0.0;
1194 pnode->Set(TO_ERASE,
false);
1195 pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=0.0;
1196 pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1199 void InterpolateOnEdge( Geometry<Node >& geom,
unsigned int step_data_size, Node::Pointer pnode)
1201 unsigned int buffer_size = pnode->GetBufferSize();
1208 double* step_data = (pnode)->SolutionStepData().Data(
step);
1211 double* node0_data = geom[0].SolutionStepData().Data(
step);
1212 double* node1_data = geom[1].SolutionStepData().Data(
step);
1215 for(
unsigned int j= 0;
j<step_data_size;
j++)
1218 step_data[
j] = 0.5*(node0_data[
j] + node1_data[
j]);
1231 pnode->FastGetSolutionStepValue(IS_LAGRANGIAN_INLET)=0.0;
1232 pnode->FastGetSolutionStepValue(IS_BOUNDARY)=1.0;
1233 pnode->FastGetSolutionStepValue(IS_FREE_SURFACE)=1.0;
1234 pnode->FastGetSolutionStepValue(IS_FLUID)=1.0;
1235 pnode->FastGetSolutionStepValue(IS_STRUCTURE)=0.0;
1238 pnode->Set(TO_ERASE,
false);
1240 if (pnode->FastGetSolutionStepValue(IS_INTERFACE)>0.99)
1241 pnode->FastGetSolutionStepValue(IS_INTERFACE)=1.0;
1243 pnode->FastGetSolutionStepValue(IS_INTERFACE)=0.0;
1251 void initialize_triangulateio( triangulateio& tr )
1253 tr.pointlist = (
REAL*) NULL;
1254 tr.pointattributelist = (
REAL*) NULL;
1255 tr.pointmarkerlist = (
int*) NULL;
1256 tr.numberofpoints = 0;
1257 tr.numberofpointattributes = 0;
1258 tr.trianglelist = (
int*) NULL;
1259 tr.triangleattributelist = (
REAL*) NULL;
1260 tr.trianglearealist = (
REAL*) NULL;
1261 tr.neighborlist = (
int*) NULL;
1262 tr.numberoftriangles = 0;
1263 tr.numberofcorners = 3;
1264 tr.numberoftriangleattributes = 0;
1265 tr.segmentlist = (
int*) NULL;
1266 tr.segmentmarkerlist = (
int*) NULL;
1267 tr.numberofsegments = 0;
1268 tr.holelist = (
REAL*) NULL;
1269 tr.numberofholes = 0;
1270 tr.regionlist = (
REAL*) NULL;
1271 tr.numberofregions = 0;
1272 tr.edgelist = (
int*) NULL;
1273 tr.edgemarkerlist = (
int*) NULL;
1274 tr.normlist = (
REAL*) NULL;
1275 tr.numberofedges = 0;
1278 void clean_triangulateio( triangulateio& tr )
1280 if(tr.pointlist != NULL) free(tr.pointlist );
1281 if(tr.pointattributelist != NULL) free(tr.pointattributelist );
1282 if(tr.pointmarkerlist != NULL) free(tr.pointmarkerlist );
1283 if(tr.trianglelist != NULL) free(tr.trianglelist );
1284 if(tr.triangleattributelist != NULL) free(tr.triangleattributelist );
1285 if(tr.trianglearealist != NULL) free(tr.trianglearealist );
1286 if(tr.neighborlist != NULL) free(tr.neighborlist );
1287 if(tr.segmentlist != NULL) free(tr.segmentlist );
1288 if(tr.segmentmarkerlist != NULL) free(tr.segmentmarkerlist );
1289 if(tr.holelist != NULL) free(tr.holelist );
1290 if(tr.regionlist != NULL) free(tr.regionlist );
1291 if(tr.edgelist != NULL) free(tr.edgelist );
1292 if(tr.edgemarkerlist != NULL) free(tr.edgemarkerlist );
1293 if(tr.normlist != NULL) free(tr.normlist );
1342 rOStream << std::endl;
PeriodicInterfaceProcess & operator=(const PeriodicInterfaceProcess &)=delete
Short class definition.
Definition: bucket.h:57
Base class for all Conditions.
Definition: condition.h:59
Geometry< NodeType >::PointsArrayType NodesArrayType
definition of nodes container type, redefined from GeometryType
Definition: condition.h:86
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
ConditionIterator ConditionsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1361
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
ConditionIterator ConditionsEnd(IndexType ThisIndex=0)
Definition: model_part.h:1371
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
void reserve(size_type dim)
Definition: pointer_vector.h:319
static double ElapsedSeconds(const std::chrono::steady_clock::time_point StartTime)
This method returns the resulting time.
Definition: timer.h:179
A generic tree data structure for spatial partitioning.
Definition: tree.h:190
Short class definition.
Definition: trigen_glass_forming.h:61
std::vector< PointType::Pointer > PointVector
Definition: trigen_glass_forming.h:70
PointVector::iterator PointIterator
Definition: trigen_glass_forming.h:71
virtual void PrintData(std::ostream &rOStream) const
Print object's data.
Definition: trigen_glass_forming.h:645
Tree< KDTreePartition< BucketType > > KdtreeType
Definition: trigen_glass_forming.h:75
void ReGenerateMesh(ModelPart &ThisModelPart, Element const &rReferenceElement, Condition const &rReferenceBoundaryCondition, EntitiesEraseProcess< Node > &node_erase, bool rem_nodes=true, bool add_nodes=true, double my_alpha=1.4, double h_factor=0.5)
Definition: trigen_glass_forming.h:105
TriGenGLASSModeler()
Default constructor.
Definition: trigen_glass_forming.h:82
virtual std::string Info() const
Turn back information as a string.
Definition: trigen_glass_forming.h:636
std::vector< double > DistanceVector
Definition: trigen_glass_forming.h:72
std::vector< double >::iterator DistanceIterator
Definition: trigen_glass_forming.h:73
KRATOS_CLASS_POINTER_DEFINITION(TriGenGLASSModeler)
Pointer definition of TriGenModeler.
virtual void PrintInfo(std::ostream &rOStream) const
Print information about this object.
Definition: trigen_glass_forming.h:642
Bucket< 3, PointType, PointVector, PointPointerType, PointIterator, DistanceIterator > BucketType
Definition: trigen_glass_forming.h:74
Node PointType
Definition: trigen_glass_forming.h:68
Node::Pointer PointPointerType
Definition: trigen_glass_forming.h:69
virtual ~TriGenGLASSModeler()
Destructor.
Definition: trigen_glass_forming.h:90
A three node 2D triangle geometry with linear shape functions.
Definition: triangle_2d_3.h:74
#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
Kratos::ModelPart ModelPart
Definition: kratos_wrapper.h:31
z
Definition: GenerateWind.py:163
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
bool CalculatePosition(Geometry< Node > &geom, const double xc, const double yc, const double zc, array_1d< double, 3 > &N)
Definition: transfer_utility.h:343
double CalculateVol(const double x0, const double y0, const double x1, const double y1, const double x2, const double y2)
Definition: transfer_utility.h:424
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
std::istream & operator>>(std::istream &rIStream, LinearMasterSlaveConstraint &rThis)
input stream function
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
void triangulate(char *, struct triangulateio *, struct triangulateio *, struct triangulateio *)
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
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
h
Definition: generate_droplet_dynamics.py:91
x2
Definition: generate_frictional_mortar_condition.py:122
x1
Definition: generate_frictional_mortar_condition.py:121
int tol
Definition: hinsberg_optimization.py:138
R
Definition: isotropic_damage_automatic_differentiation.py:172
float radius
Definition: mesh_to_mdpa_converter.py:18
def Interpolate(variable, entity, sf_values, historical_value)
Definition: point_output_process.py:231
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
namespace
Definition: array_1d.h:793
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31
float factor
for node in (self.combined_model_part).Nodes: pold = node.GetSolutionStepValue(PRESSURE,...
Definition: ulf_PGLASS.py:254