22 plucker_coordinates_abscissa_ =
nullptr;
23 plucker_coordinates_triangle_.resize(3,
nullptr);
24 plucker_products_.resize(3,
nullptr);
35 plucker_coordinates_abscissa_ =
new Plucker(*abscissa.
node(0),
36 *abscissa.
node(1),
true);
37 scale_line_=plucker_coordinates_abscissa_->scale();
39 plucker_coordinates_triangle_.resize(3);
40 plucker_products_.resize(3);
41 scale_triangle_=std::numeric_limits<double>::max();
42 for(
unsigned int side = 0; side < 3; side++){
46 scale_triangle_ = std::min( scale_triangle_, plucker_coordinates_triangle_[side]->scale());
49 plucker_products_[side] =
new double((*plucker_coordinates_abscissa_)*(*plucker_coordinates_triangle_[side]));
55 if(plucker_coordinates_abscissa_ !=
nullptr)
56 delete plucker_coordinates_abscissa_;
58 for(
unsigned int side = 0; side < RefElement<2>::n_sides; side++){
59 if(plucker_products_[side] !=
nullptr)
60 delete plucker_products_[side];
61 if(plucker_coordinates_triangle_[side] !=
nullptr)
62 delete plucker_coordinates_triangle_[side];
69 for(
unsigned int side = 0; side < RefElement<2>::n_sides; side++)
71 plucker_products_[side] =
nullptr;
72 plucker_coordinates_triangle_[side] =
nullptr;
74 plucker_coordinates_abscissa_ =
nullptr;
80 plucker_coordinates_abscissa_->compute();
81 scale_line_=plucker_coordinates_abscissa_->scale();
86 scale_triangle_=std::numeric_limits<double>::max();
88 for(
unsigned int side = 0; side < RefElement<2>::n_sides; side++){
89 plucker_coordinates_triangle_[side]->compute();
90 scale_triangle_ = std::min( scale_triangle_, plucker_coordinates_triangle_[side]->scale());
92 ASSERT(plucker_products_[side]).error(
"Undefined plucker product.");
94 *plucker_products_[side] = (*plucker_coordinates_abscissa_)*(*plucker_coordinates_triangle_[side]);
107 int sign = (t < tol ? -2 : (t > 1+tol ? +2 : 0) );
108 if (std::abs(t) <= tol){
112 else if (std::abs(1-t) <= tol){
141 (local[0])(local[1])(local[2]);
145 arma::vec3 local_triangle({local[2],local[1],local[0]});
151 arma::vec3 u = plucker_coordinates_abscissa_->get_u_vector();
156 if(fabs((
double)u[1]) > fabs(max)){ max = u[1]; i = 1;}
157 if(fabs((
double)u[2]) > fabs(max)){ max = u[2]; i = 2;}
160 double isect_coord_i =
166 double t = (-plucker_coordinates_abscissa_->point(0)[i] + isect_coord_i)/max;
168 arma::vec2 local_abscissa = {1-t, t};
187 compute_plucker_products();
196 signed_plucker_product(1),
197 signed_plucker_product(2)};
198 double w_sum = w[0] + w[1] + w[2];
200 unsigned int n_positive = 0;
201 unsigned int n_negative = 0;
202 unsigned int zero_idx_sum =0;
207 double scaled_epsilon =
geometry_epsilon*scale_line_*scale_triangle_*scale_triangle_;
208 if(std::fabs(w_sum) > scaled_epsilon) {
210 for (
unsigned int i=0; i < 3; i++) {
225 for (
unsigned int i=0; i < 3; i++) {
227 if (w[i] > scaled_epsilon || w[i] < -scaled_epsilon) n_negative++;
239 if (n_positive > 0) {
241 compute_plucker(IP, w);
243 unsigned int non_zero_idx=0;
244 if (n_positive == 2) {
248 non_zero_idx = (zero_idx_sum + 1) % 3;
250 else if (n_positive == 1) {
254 non_zero_idx = 3-zero_idx_sum;
282 int sign = check_abscissa_topology(IP);
284 if(std::abs(sign) > 1)
return 0;
292 return compute_final_in_plane(IP12s);
320 arma::vec3 A = plucker_coordinates_abscissa_->point(0);
322 arma::vec3 U = plucker_coordinates_abscissa_->get_u_vector();
326 arma::vec3 V = plucker_coordinates_triangle_[side]->get_u_vector();
344 unsigned int max_index = 0;
345 double maximum = fabs(Det[0]);
346 if(fabs((
double)Det[1]) > maximum){
347 maximum = fabs(Det[1]);
350 if(fabs((
double)Det[2]) > maximum){
351 maximum = fabs(Det[2]);
365 unsigned int i = (max_index+1)%3,
368 double DetX = -K[i]*V[j] + K[j]*V[i];
369 double DetY = -K[i]*U[j] + K[j]*U[i];
371 double s = DetX/Det[max_index];
372 double t = DetY/Det[max_index];
383 if(t >= -tol && t <= 1+tol){
392 arma::vec2 local_abscissa({1-s, s});
413 for(
unsigned int i = 0; i < 3;i++){
415 if(IP12s.size() == 2)
break;
418 if (compute_degenerate(i,IP))
420 uint s = check_abscissa_topology(IP);
429 if(IP12s.size() == 0)
return IP12s.size();
432 if(IP12s.size() == 1){
433 if(std::abs(sign[0]) > 1) IP12s.pop_back();
446 if (sign[0] == sign[1] && std::abs(sign[0]) > 1) {
452 if(IP12s[0].local_bcoords_A()[1] > IP12s[1].local_bcoords_A()[1]){
458 const int ip_sign[] = {-2, +2};
459 for(
unsigned int ip=0; ip<2; ip++) {
461 if (sign[ip] == ip_sign[ip]) {
463 correct_triangle_ip_topology(
double(ip), ip, IP12s);
464 IP12s[ip].set_topology_A(ip, 0);
469 if(IP12s[0].topology_equal(IP12s[1])){
481 ips[0].local_bcoords_B(), ips[1].local_bcoords_B(),
482 ips[0].local_bcoords_A()[1], ips[1].local_bcoords_A()[1], t);
483 arma::vec2 local_abscissa({1 - t, t});
484 ips[ip].set_coordinates(local_abscissa, local_tria);
489 std::pair<unsigned int, unsigned int> zeros =
502 default: ips[ip].set_topology_B(0,2);
513 os <<
"\tPluckerCoordinates Abscissa[0]";
514 if(plucker_coordinates_abscissa_ ==
nullptr){
515 os <<
"NULL" << endl;
517 os << plucker_coordinates_abscissa_;
519 for(
unsigned int i = 0; i < 3;i++){
520 os <<
"\tPluckerCoordinates Triangle[" << i <<
"]";
521 if(plucker_coordinates_triangle_[i] ==
nullptr){
522 os <<
"NULL" << endl;
524 os << plucker_coordinates_triangle_[i];
549 for(
unsigned int side = 0; side < RefElement<2>::n_sides; side++){
558 for(
unsigned int p = 0; p < 3*RefElement<2>::n_sides; p++){
566 for(
unsigned int side = 0; side < 2*RefElement<2>::n_sides; side++)
567 CI12[side].clear_all();
570 for(
unsigned int side = 0; side < 2*RefElement<2>::n_sides; side++){
571 if(plucker_coordinates_[side] !=
nullptr)
572 delete plucker_coordinates_[side];
575 for(
unsigned int p = 0; p < 3*RefElement<2>::n_sides; p++){
576 if(plucker_products_[p] !=
nullptr)
577 delete plucker_products_[p];
584 for(
unsigned int side = 0; side < 2*RefElement<2>::n_sides; side++)
586 plucker_coordinates_[side] =
nullptr;
588 for(
unsigned int p = 0; p < 3*RefElement<2>::n_sides; p++){
589 plucker_products_[p] =
nullptr;
596 for(
unsigned int i = 0; i < RefElement<2>::n_sides; i++){
598 for(
unsigned int j = 0; j < RefElement<2>::n_sides; j++)
599 CI12[i].set_pc_triangle(plucker_coordinates_[3+j], j);
601 CI12[i].set_pc_abscissa(plucker_coordinates_[i]);
604 for(
unsigned int j = 0; j < RefElement<2>::n_sides; j++)
610 for(
unsigned int j = 0; j < RefElement<2>::n_sides; j++)
613 CI12[i].set_plucker_product(get_plucker_product(i,j),j);
629 unsigned int ip_coutner = 0;
632 for(
unsigned int i = 0; i < 2*RefElement<2>::n_sides && ip_coutner < 2; i++){
633 if(!CI12[i].is_computed())
636 if(CI12[i].compute_final(IP12s) == 0)
continue;
638 unsigned int triangle_side = i%3;
642 IPAux22 IP22(IP.switch_objects(), triangle_side);
648 if( IP22.
dim_A() == 0 )
655 for(
unsigned int s=0; s < RefElement<2>::n_sides_per_node; s++)
658 if( IP22.
dim_B() == 0 )
662 for(
unsigned int s=0; s < RefElement<2>::n_sides_per_node; s++)
665 else if( IP22.
dim_B() == 1 )
669 CI12[3 + IP22.
idx_B()].set_computed();
672 else if( IP22.
dim_B() == 0 )
680 for(
unsigned int s=0; s < RefElement<2>::n_sides_per_node; s++)
685 IP22s.push_back(IP22);
723 for(
unsigned int i = 0; i < RefElement<2>::n_lines; i++){
724 os <<
"\tPluckerCoordinates Triangle A[" << i <<
"]";
725 if(plucker_coordinates_[i] ==
nullptr){
726 os <<
"NULL" << endl;
728 os << plucker_coordinates_[i];
731 for(
unsigned int i = 0; i < RefElement<2>::n_lines; i++){
732 os <<
"\tPluckerCoordinates Triangle B[" << i <<
"]";
734 os <<
"NULL" << endl;
736 os << plucker_coordinates_[RefElement<2>::n_lines + i];
742 os <<
"ComputeIntersection<2,2> Plucker Coordinates Tree:" << endl;
743 print_plucker_coordinates(os);
744 for(
unsigned int i = 0; i < 6;i++){
745 os <<
"ComputeIntersection<1,2>["<< i <<
"] Plucker Coordinates:" << endl;
746 CI12[i].print_plucker_coordinates(os);
756 plucker_coordinates_abscissa_ =
nullptr;
757 plucker_coordinates_tetrahedron.resize(6,
nullptr);
758 plucker_products_.resize(6,
nullptr);
767 "Tetrahedron element (%d) has wrong numbering or is degenerated (negative Jacobian).");
770 plucker_coordinates_abscissa_ =
new Plucker(*abscissa.
node(0), *abscissa.
node(1));
771 plucker_coordinates_tetrahedron.resize(6);
772 plucker_products_.resize(6);
774 for(
unsigned int line = 0; line < RefElement<3>::n_lines; line++){
785 for(
unsigned int side = 0; side < RefElement<3>::n_sides; side++)
786 CI12[side].clear_all();
789 if(plucker_coordinates_abscissa_ !=
nullptr)
790 delete plucker_coordinates_abscissa_;
792 for(
unsigned int line = 0; line < RefElement<3>::n_lines; line++){
793 if(plucker_products_[line] !=
nullptr)
794 delete plucker_products_[line];
795 if(plucker_coordinates_tetrahedron[line] !=
nullptr)
796 delete plucker_coordinates_tetrahedron[line];
803 for(
unsigned int side = 0; side < RefElement<3>::n_lines; side++)
805 plucker_products_[side] =
nullptr;
806 plucker_coordinates_tetrahedron[side] =
nullptr;
808 plucker_coordinates_abscissa_ =
nullptr;
814 for(
unsigned int side = 0; side < RefElement<3>::n_sides; side++){
815 for(
unsigned int line = 0; line < RefElement<3>::n_lines_per_side; line++){
816 CI12[side].set_pc_triangle(plucker_coordinates_tetrahedron[
819 CI12[side].set_plucker_product(
823 CI12[side].set_pc_abscissa(plucker_coordinates_abscissa_);
841 for(
unsigned int face = 0; face < RefElement<3>::n_sides && IP13s.size() < 2; face++){
843 if (CI12[face].is_computed())
continue;
851 position_on_abscissa.push_back(CI12[face].check_abscissa_topology(IP12));
853 IPAux IP13(IP12, face);
855 IP13s.push_back(IP13);
858 if(IP13.
dim_B() == 0)
862 CI12[node_face].set_computed();
864 else if(IP13.
dim_B() == 1)
867 CI12[edge_face].set_computed();
871 if (IP13s.size() == 0)
return IP13s.size();
874 if(IP13s.size() == 1){
875 if(std::abs(position_on_abscissa[0]) > 1) IP13s.pop_back();
883 if (position_on_abscissa[0] == position_on_abscissa[1] && std::abs(position_on_abscissa[0]) > 1) {
889 if(IP13s[0].local_bcoords_A()[1] > IP13s[1].local_bcoords_A()[1]){
891 std::swap(position_on_abscissa[0], position_on_abscissa[1]);
895 const int ip_sign[] = {-2, +2};
896 for(
unsigned int ip=0; ip<2; ip++) {
898 if (position_on_abscissa[ip] == ip_sign[ip]) {
899 position_on_abscissa[ip] /=2;
900 correct_tetrahedron_ip_topology(
double(ip), ip, IP13s);
901 IP13s[ip].set_topology_A(ip, 0);
906 if(IP13s[0].topology_equal(IP13s[1])){
917 ips[0].local_bcoords_B(), ips[1].local_bcoords_B(),
918 ips[0].local_bcoords_A()[1], ips[1].local_bcoords_A()[1], t);
919 arma::vec2 local_abscissa({1 - t, t});
920 ips[ip].set_coordinates(local_abscissa, local_tetra);
925 unsigned int zeros = 0;
926 unsigned int n_zeros = 0;
927 for(
char i=0; i < 4; i++){
930 zeros = zeros | (1 << i);
944 default: ips[ip].set_topology_B(0,3);
957 os <<
"\tPluckerCoordinates Abscissa[0]";
958 if(plucker_coordinates_abscissa_ ==
nullptr){
959 os <<
"NULL" << endl;
961 os << plucker_coordinates_abscissa_;
964 for(
unsigned int i = 0; i < 6;i++){
965 os <<
"\tPluckerCoordinates Tetrahedron[" << i <<
"]";
966 if(plucker_coordinates_tetrahedron[i] ==
nullptr){
967 os <<
"NULL" << endl;
969 os << plucker_coordinates_tetrahedron[i];
975 os <<
"ComputeIntersection<1,3> Plucker Coordinates Tree:" << endl;
976 print_plucker_coordinates(os);
977 for(
unsigned int i = 0; i < 4;i++){
978 os <<
"ComputeIntersection<1,2>["<< i <<
"] Plucker Coordinates:" << endl;
979 CI12[i].print_plucker_coordinates(os);
990 s3_dim_starts({0, 4, 10, 14}),
991 s2_dim_starts({15, 18, 21}),
992 object_next(22, no_idx),
993 on_faces(_on_faces())
996 plucker_coordinates_triangle_.resize(3,
nullptr);
997 plucker_coordinates_tetrahedron.resize(6,
nullptr);
998 plucker_products_.resize(3*6,
nullptr);
1008 "Tetrahedron element (%d) has wrong numbering or is degenerated (negative Jacobian).");
1010 S3_inverted = tetrahedron.
inverted();
1012 plucker_coordinates_triangle_.resize(3);
1013 plucker_coordinates_tetrahedron.resize(6);
1016 for(
unsigned int i = 0; i < RefElement<3>::n_lines; i++){
1021 for(
unsigned int i = 0; i < RefElement<2>::n_lines;i++){
1030 plucker_products_.resize(np,
nullptr);
1031 for(
unsigned int line = 0; line < np; line++){
1040 for(
unsigned int triangle_side = 0; triangle_side < RefElement<2>::n_sides; triangle_side++)
1041 CI13[triangle_side].clear_all();
1043 for(
unsigned int line = 0; line < RefElement<3>::n_lines; line++)
1044 CI12[line].clear_all();
1048 for(
unsigned int line = 0; line < np; line++){
1049 if(plucker_products_[line] !=
nullptr)
1050 delete plucker_products_[line];
1053 for(
unsigned int i = 0; i < RefElement<3>::n_lines;i++){
1054 if(plucker_coordinates_tetrahedron[i] !=
nullptr)
1055 delete plucker_coordinates_tetrahedron[i];
1057 for(
unsigned int i = 0; i < RefElement<2>::n_sides;i++){
1058 if(plucker_coordinates_triangle_[i] !=
nullptr)
1059 delete plucker_coordinates_triangle_[i];
1068 for(
unsigned int triangle_side = 0; triangle_side < RefElement<2>::n_sides; triangle_side++){
1069 for(
unsigned int line = 0; line < RefElement<3>::n_lines; line++){
1070 CI13[triangle_side].set_plucker_product(
1073 CI12[line].set_plucker_product(
1077 CI13[triangle_side].set_pc_tetrahedron(plucker_coordinates_tetrahedron[line],line);
1078 CI12[line].set_pc_triangle(plucker_coordinates_triangle_[triangle_side],triangle_side);
1080 CI13[triangle_side].set_pc_abscissa(plucker_coordinates_triangle_[triangle_side]);
1081 CI13[triangle_side].init();
1085 for(
unsigned int line = 0; line < RefElement<3>::n_lines; line++)
1086 CI12[line].set_pc_abscissa(plucker_coordinates_tetrahedron[line]);
1094 unsigned int ip = object_next[i_obj];
1095 if (ip == no_idx)
return false;
1097 return IP_next[ip] == i_obj;
1105 if (have_backlink(obj_after_ip)) {
1110 unsigned int ip = object_next[obj_after_ip];
1111 ASSERT( ! have_backlink(obj_after_ip) )
1112 (intersection_->component_ele_idx())
1113 (intersection_->bulk_ele_idx())
1118 (object_next[obj_after_ip]);
1121 object_next[obj_before_ip] = ip_idx;
1122 IP_next.push_back( obj_after_ip);
1123 if (object_next[obj_after_ip] == no_idx) {
1124 object_next[obj_after_ip] = ip_idx;
1132 intersection_= &intersection;
1139 std::fill(object_next.begin(), object_next.end(), no_idx);
1142 std::array<bool, 6> edge_touch={
false,
false,
false,
false,
false,
false};
1143 unsigned int object_before_ip, object_after_ip;
1162 for(
unsigned int _i_side = 0; _i_side < RefElement<2>::n_lines; _i_side++) {
1163 unsigned int i_side = cycle_sides[_i_side];
1165 CI13[ i_side ].compute(IP13s);
1166 ASSERT(IP13s.size() < 3);
1167 if (IP13s.size() == 0)
continue;
1168 for(
unsigned int _ip=0; _ip < IP13s.size(); _ip++) {
1172 unsigned int ip = (side_cycle_orientation[_i_side] + _ip) % IP13s.size();
1183 unsigned int tetra_object = s3_dim_starts[IP23.
dim_B()] + IP23.
idx_B();
1184 unsigned int side_object = s2_dim_starts[1] + i_side;
1186 object_before_ip = tetra_object;
1187 object_after_ip = side_object;
1195 object_before_ip = s2_dim_starts[0] + IP23.
idx_A();
1199 if (IP13s.size() == 1 ) {
1200 if (IP23.
dim_B() == 0) {
1203 if (IP23.
dim_A() == 0) {
1204 object_before_ip = tetra_object;
1205 object_after_ip = s2_dim_starts[0]+IP23.
idx_A();
1215 edge_touch[IP23.
idx_B()]=
true;
1216 std::swap(object_before_ip, object_after_ip);
1220 unsigned int ip_idx = add_ip(IP23);
1223 ASSERT_EQ(IP23_list.size(), IP_next.size()+1);
1224 set_links(object_before_ip, ip_idx, object_after_ip);
1244 for(
unsigned int tetra_edge = 0; tetra_edge < 6; tetra_edge++) {
1246 CI12[tetra_edge].compute(IP12);
1250 IP12s_.push_back(IP12);
1254 for(
unsigned int tetra_edge = 0; tetra_edge < 6;tetra_edge++) {
1255 if (! processed_edge[tetra_edge]) {
1257 IPAux12 &IP12 = IP12s_[tetra_edge];
1260 int sign = CI12[tetra_edge].check_abscissa_topology(IP12);
1262 if(std::abs(sign) > 1)
continue;
1264 IPAux23 IP23(IP12s_[tetra_edge].switch_objects(), tetra_edge);
1270 unsigned int ip_idx = add_ip(IP23);
1273 if ( edge_dim == 0) {
1274 face_pair = vertex_faces(i_edge);
1277 processed_edge[ie] = 1;
1280 face_pair = edge_faces(i_edge);
1286 unsigned int s3_object = s3_dim_starts[edge_dim] + i_edge;
1289 if (IP23.
dim_A() < 2
1290 && (! edge_touch[i_edge])
1291 && object_next[s3_object] != no_idx) {
1293 if ( have_backlink(s3_object) ) {
1294 set_links(s3_object, ip_idx, face_pair[1]);
1296 set_links(face_pair[0], ip_idx, s3_object);
1301 set_links(face_pair[0], ip_idx, face_pair[1]);
1303 if ( have_backlink(s3_object) ) {
1304 object_next[s3_object]=ip_idx;
1314 if (IP23_list.size() == 0)
return;
1320 for(
auto obj : IP_next) {
1322 unsigned int ip = object_next[obj];
1323 if (ip < IP_next.size()) have_predecessor[ip]=1;
1325 unsigned int ip_init=0;
1326 for(
unsigned int i=0; i< IP23_list.size(); i++)
if (! have_predecessor[i]) ip_init=i;
1329 arma::uvec::fixed<RefElement<3>::n_sides> ips_face_counter;
1330 ips_face_counter.zeros();
1333 unsigned int ip=ip_init;
1334 ASSERT_EQ(IP_next.size(), IP23_list.size());
1335 intersection.
points().push_back(IP23_list[ip]);
1340 unsigned int object = IP_next[ip];
1343 ip = object_next[object];
1344 object_next[object]=no_idx;
1345 if ((ip == no_idx))
break;
1348 if ( ! IP23_list[ip].topology_equal(intersection.
points().back()) ) {
1351 intersection.
points().push_back(IP);
1353 ips_face_counter += on_faces[IP.
dim_B()][IP.
idx_B()];
1357 if (intersection.
points().size() == 1)
return;
1359 if (IP23_list[ip_init].topology_equal(intersection.
points().back()) )
1360 intersection.
points().pop_back();
1362 if (intersection.
points().size() > 2){
1363 for(
uint i=0; i< ips_face_counter.n_elem; i++){
1364 if(ips_face_counter(i) >= intersection.
points().size()){
1376 unsigned int ip_ori = (
unsigned int)(IP12s_[i_edge].result());
1383 return { s3_dim_starts[2] + line_faces[1-ip_ori], s3_dim_starts[2] + line_faces[ip_ori] };
1390 std::array<unsigned int, 3> n_ori, sum_idx;
1393 for(
unsigned int ie=0; ie <3; ie++) {
1394 unsigned int edge_ip_ori = (
unsigned int)(IP12s_[ vtx_edges[ie]].result());
1397 edge_ip_ori = (edge_ip_ori +1)%2;
1399 if (edge_ip_ori == 3) edge_ip_ori=2;
1400 n_ori[edge_ip_ori]++;
1401 sum_idx[edge_ip_ori]+=ie;
1408 if ( n_degen == 2 ) {
1412 unsigned int i_edge = 3 - sum_degen;
1413 FacePair pair = edge_faces(vtx_edges[i_edge]);
1416 if (pair[0] == s3_dim_starts[2] + vtx_faces[(i_edge+1)%3])
1417 return { s3_dim_starts[1] + (i_edge+2)%3, s3_dim_starts[1] + (i_edge+1)%3 };
1419 return { s3_dim_starts[1] + (i_edge+1)%3, s3_dim_starts[1] + (i_edge+2)%3 };
1421 }
else if (n_degen == 1) {
1423 unsigned int i_edge = sum_degen;
1424 ASSERT( n_positive + n_negative == 2);
1425 if ( n_positive == 1) {
1428 FacePair pair = edge_faces(vtx_edges[(i_edge+1)%3]);
1433 if (pair[0] == s3_dim_starts[2] + face)
1434 return { s3_dim_starts[2] + face, s3_dim_starts[1] + vtx_edges[i_edge]};
1436 return { s3_dim_starts[1] + vtx_edges[i_edge], s3_dim_starts[2] + face };
1441 ASSERT(n_positive == 0 || n_positive== 2);
1442 return { s3_dim_starts[0]+i_vertex, s3_dim_starts[1] + vtx_edges[i_edge]};
1448 ASSERT( n_positive + n_negative == 3);
1450 if (n_positive == 1) {
1452 return edge_faces(vtx_edges[i_edge]);
1453 }
else if (n_negative == 1) {
1455 return edge_faces(vtx_edges[i_edge]);
1458 ASSERT( n_positive == 0 || n_positive == 3);
1459 return { s3_dim_starts[0]+i_vertex, s3_dim_starts[0]+i_vertex};
1470 arma::uvec::fixed<RefElement<3>::n_sides> v; v.zeros();
1473 for(
uint i=0; i<RefElement<3>::n_nodes; i++) {
1475 for(
uint j=0; j<RefElement<3>::n_sides_per_node; j++)
1476 on_faces[0][i](faces[j]) = 1;
1480 for(
uint i=0; i<RefElement<3>::n_lines; i++) {
1482 for(
uint j=0; j<RefElement<3>::n_sides_per_line; j++)
1483 on_faces[1][i](faces[j]) = 1;
1487 for(
uint i=0; i<RefElement<3>::n_sides; i++) {
1488 on_faces[2][i](i) = 1;
1503 for(
unsigned int i = 0; i < 3;i++){
1504 os <<
"\tPluckerCoordinates Triangle[" << i <<
"]";
1505 if(plucker_coordinates_triangle_[i] ==
nullptr){
1506 os <<
"NULL" << endl;
1508 os << plucker_coordinates_triangle_[i];
1511 for(
unsigned int i = 0; i < 6;i++){
1512 os <<
"\tPluckerCoordinates Tetrahedron[" << i <<
"]";
1513 if(plucker_coordinates_tetrahedron[i] ==
nullptr){
1514 os <<
"NULL" << endl;
1516 os << plucker_coordinates_tetrahedron[i];
1522 os <<
"ComputeIntersection<2,3> Plucker Coordinates Tree:" << endl;
1523 print_plucker_coordinates(os);
1524 for(
unsigned int i = 0; i < 6;i++){
1525 os <<
"ComputeIntersection<1,2>["<< i <<
"] Plucker Coordinates:" << endl;
1526 CI12[i].print_plucker_coordinates(os);
1528 for(
unsigned int i = 0; i < 3;i++){
1529 CI13[i].print_plucker_coordinates_tree(os);
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
#define ASSERT_PERMANENT_EQ(a, b)
Definition of comparative assert macro (EQual)
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
std::array< uint, 2 > FacePair
unsigned int input_id() const
Return the element ID in the input mesh. Should be only used for special output.
NodeAccessor< 3 > node(unsigned int ni) const
double jacobian_S3() const
Internal auxiliary class representing intersection object of simplex<dimA> and simplex<dimB>.
void set_ips_in_face(unsigned int face_idx)
unsigned int size() const
Returns number of intersection points.
std::vector< IntersectionPointAux< dimA, dimB > > i_points_
Vector of internal intersection points.
std::vector< IntersectionPointAux< dimA, dimB > > & points()
Returns intersection points by a reference.
Internal auxiliary class represents an intersection point of simplex<N> and simplex<M>.
unsigned int idx_B() const
Returns the index of Simplex<M>.
bool topology_equal(const IntersectionPointAux< N, M > &other) const
Returns true, if other intersection point has the same topology.
void set_coordinates(const arma::vec::fixed< N+1 > &lcA, const arma::vec::fixed< M+1 > &lcB)
Setter for coordinates.
IntersectionResult result() const
Result: 0 - negative sign, 1 - positive sign, 2 - degenerate (zero for all sides),...
void set_result(IntersectionResult result)
Setter orientation flag.
unsigned int dim_A() const
Returns dimension of object A.
unsigned int idx_A() const
Returns the index of Simplex<N>.
void set_topology_B(unsigned int idx, unsigned int dim_B)
Sets the topology of object B in Simplex<M>.
void set_topology_A(unsigned int idx, unsigned int dim_A)
Sets the topology of object A in Simplex<N>.
unsigned int dim_B() const
Returns dimension of object B.
const arma::vec::fixed< N+1 > & local_bcoords_A() const
Returns barycentric coordinates in the Simplex<N>.
IntersectionPointAux< M, N > switch_objects() const
Switches the object A and B.
Plucker coordinates representing line given by points A,B.
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
static unsigned int normal_orientation(unsigned int sid)
static BaryPoint line_barycentric_interpolation(BaryPoint first_coords, BaryPoint second_coords, double first_theta, double second_theta, double theta)
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon() *2)
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact(TInteraction< OutDim, InDim > interaction, bool inv=false)
Fundamental simplicial intersections.
static const double plucker_empty
Auxiliary value for Plucker product. If equal this value, it is supposed not to be set yet.
Internal class representing intersection object.
Internal class representing intersection point.
static const double geometry_epsilon
void swap(nlohmann::json &j1, nlohmann::json &j2) noexcept(is_nothrow_move_constructible< nlohmann::json >::value and is_nothrow_move_assignable< nlohmann::json >::value)
exchanges the values of two JSON objects
Plucker coordinates class.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
std::array< unsigned int, Size > IdxVector