10#pragma mark * local ( static ) function prototypes *
12static double CurrentTime(
void)
27#pragma mark * exported function implementations *
68 for(
int i=0,eN=0;i<nx-1;i++)
75 const double hx=Lx/(nx-1.0);
77 memset(mesh.
nodeArray,0,nx*3*
sizeof(
double));
99 for(
int i=0,eN=0;i<ny-1;i++)
100 for(
int j=0;j<nx-1;j++)
103 n0 =(j+0) + (i+0)*nx,
104 n1 =(j+1) + (i+0)*nx,
105 n2 =(j+0) + (i+1)*nx,
106 n3 =(j+1) + (i+1)*nx;
107 if (triangleFlag == 2)
113 else if (triangleFlag == 1)
115 if (i%2 + j%2 == 0 || i%2 + j%2 == 2)
144 double x_0,y_0,x_1,y_1,epsilon=1.0e-8;
152 if (y_0 <= epsilon && y_1 <= epsilon)
154 else if (y_0 >= Ly - epsilon && y_1 >= Ly - epsilon)
156 else if (x_0 <= epsilon && x_1 <= epsilon)
158 else if (x_0 >= Lx - epsilon && x_1 >= Lx - epsilon)
168 const double hx=Lx/(nx-1.0),hy=Ly/(ny-1.0);
177 for(
int i=0;i<ny;i++)
178 for(
int j=0;j<nx;j++)
199 int &n0(nodes[0]),&n1(nodes[1]),&n2(nodes[2]),&n3(nodes[3]);
201 t[0][0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
202 t[0][1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
203 t[0][2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
205 t[1][0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
206 t[1][1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
207 t[1][2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
209 t[2][0] = nodeArray[n3*3+0] - nodeArray[n0*3+0];
210 t[2][1] = nodeArray[n3*3+1] - nodeArray[n0*3+1];
211 t[2][2] = nodeArray[n3*3+2] - nodeArray[n0*3+2];
213 double det = t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) -
214 t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) +
215 t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]);
245 for(
int i=0,eN=0;i<nz-1;i++)
246 for(
int j=0;j<ny-1;j++)
247 for(
int k=0;k<nx-1;k++)
250 n0 = (k+0) + (j+0)*nx + (i+0)*nxy,
251 n1 = (k+1) + (j+0)*nx + (i+0)*nxy,
252 n2 = (k+0) + (j+1)*nx + (i+0)*nxy,
253 n3 = (k+1) + (j+1)*nx + (i+0)*nxy,
254 n4 = (k+0) + (j+0)*nx + (i+1)*nxy,
255 n5 = (k+1) + (j+0)*nx + (i+1)*nxy,
256 n6 = (k+0) + (j+1)*nx + (i+1)*nxy,
257 n7 = (k+1) + (j+1)*nx + (i+1)*nxy;
272 int ebN,nN_0,nN_1,nN_2;
294 if (z_0 <= epsilon && z_1 <= epsilon && z_2 <= epsilon)
296 else if (z_0 >= Lz - epsilon && z_1 >= Lz - epsilon && z_2 >= Lz - epsilon)
298 else if (y_0 <= epsilon && y_1 <= epsilon && y_2 <= epsilon)
300 else if (y_0 >= Ly - epsilon && y_1 >= Ly - epsilon && y_2 >= Ly - epsilon)
302 else if (x_0 <= epsilon && x_1 <= epsilon && x_2 <= epsilon)
304 else if (x_0 >= Lx - epsilon && x_1 >= Lx - epsilon && x_2 >= Lx - epsilon)
316 int ebN,nN_0,nN_1,nN_2, nN_3;
344 if (z_0 <= epsilon && z_1 <= epsilon && z_2 <= epsilon && z_3 <= epsilon)
346 else if (z_0 >= Lz - epsilon && z_1 >= Lz - epsilon && z_2 >= Lz - epsilon && z_3 >= Lz - epsilon)
348 else if (y_0 <= epsilon && y_1 <= epsilon && y_2 <= epsilon && y_3 <= epsilon)
350 else if (y_0 >= Ly - epsilon && y_1 >= Ly - epsilon && y_2 >= Ly - epsilon && y_3 >= Ly - epsilon)
352 else if (x_0 <= epsilon && x_1 <= epsilon && x_2 <= epsilon && x_3 <= epsilon)
354 else if (x_0 >= Lx - epsilon && x_1 >= Lx - epsilon && x_2 >= Lx - epsilon && x_3 >= Lx - epsilon)
377 const double hx=Lx/(nx-1.0),
393 for(
int i=0;i<nz;i++)
394 for(
int j=0;j<ny;j++)
395 for(
int k=0;k<nx;k++)
397 nN = k + j*nx + i*nxy;
425 const double hx=Lx/(nx-1.0),
436 for(
int j=0;j<ny;j++)
437 for(
int k=0;k<nx;k++)
478 std::cout<<
"mesh elements"<<std::endl;
496 for(
int i=0;i<nz-1;i++)
497 for(
int j=0;j<ny-1;j++)
498 for(
int k=0;k<nx-1;k++)
501 (k+0) + (j+0)*nx + (i+0)*nxy,
502 (k+1) + (j+0)*nx + (i+0)*nxy,
503 (k+1) + (j+1)*nx + (i+0)*nxy,
504 (k+0) + (j+1)*nx + (i+0)*nxy,
505 (k+0) + (j+0)*nx + (i+1)*nxy,
506 (k+1) + (j+0)*nx + (i+1)*nxy,
507 (k+1) + (j+1)*nx + (i+1)*nxy,
508 (k+0) + (j+1)*nx + (i+1)*nxy);
529 for(
int j=0;j<ny-1;j++)
530 for(
int k=0;k<nx-1;k++)
568 for(
int i=0;i<nz-px;i++)
569 for(
int j=0;j<ny-py;j++)
570 for(
int k=0;k<nx-pz;k++)
577 for(
int ii=0;ii<px+1;ii++)
578 for(
int jj=0;jj<py+1;jj++)
579 for(
int kk=0;kk<pz+1;kk++)
588 mesh.
U_KNOT =
new double[nx+px+1];
589 mesh.
V_KNOT =
new double[ny+py+1];
590 mesh.
W_KNOT =
new double[nz+pz+1];
592 for(
int i=0;i<px+1;i++)
594 for(
int i=px+1;i<nx;i++)
595 mesh.
U_KNOT[i] =
double(i-px-1);
596 for(
int i=nx;i<nx+px+1;i++)
597 mesh.
U_KNOT[i] =
double(nx);
599 for(
int i=0;i<py+1;i++)
601 for(
int i=py+1;i<ny;i++)
602 mesh.
V_KNOT[i] =
double(i-py-1);
603 for(
int i=ny;i<ny+py+1;i++)
604 mesh.
V_KNOT[i] =
double(ny);
606 for(
int i=0;i<pz+1;i++)
608 for(
int i=pz+1;i<pz;i++)
609 mesh.
W_KNOT[i] =
double(i-pz-1);
610 for(
int i=nz;i<nz+pz+1;i++)
611 mesh.
W_KNOT[i] =
double(nz);
628 logEvent(
"Constructing element boundary map",6);
633 for(
int ebN=0;ebN<2;ebN++)
638 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
640 elementBoundaryElements[ebt].right=eN;
641 elementBoundaryElements[ebt].right_ebN_element=ebN;
645 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
648 stop = CurrentTime();
654 start = CurrentTime();
655 set<int> interiorElementBoundaries,exteriorElementBoundaries;
662 stop = CurrentTime();
666 start = CurrentTime();
669 eb != elementBoundaryElements.end();
678 if(eb->second.right != -1)
680 interiorElementBoundaries.insert(ebN);
684 exteriorElementBoundaries.insert(ebN);
687 if (eb->second.right != -1)
697 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
699 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
721 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
726 stop = CurrentTime();
741 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
791 for(
int ebN=0;ebN<3;ebN++)
797 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
799 elementBoundaryElements[ebt].right=eN;
800 elementBoundaryElements[ebt].right_ebN_element=ebN;
804 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
807 stop = CurrentTime();
813 start = CurrentTime();
814 set<int> interiorElementBoundaries,exteriorElementBoundaries;
820 stop = CurrentTime();
824 start = CurrentTime();
827 eb != elementBoundaryElements.end();
837 if(eb->second.right != -1)
839 interiorElementBoundaries.insert(ebN);
843 exteriorElementBoundaries.insert(ebN);
845 if (eb->second.right != -1)
855 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
857 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
859 set<NodeTuple<2> > edges;
874 for (set<
NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
891 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
893 stop = CurrentTime();
911 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
960 for(
int ebN=0;ebN<4;ebN++)
966 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
968 elementBoundaryElements[ebt].right=eN;
969 elementBoundaryElements[ebt].right_ebN_element=ebN;
973 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
976 stop = CurrentTime();
982 start = CurrentTime();
983 set<int> interiorElementBoundaries,exteriorElementBoundaries;
989 stop = CurrentTime();
993 start = CurrentTime();
996 eb != elementBoundaryElements.end();
1006 if(eb->second.right != -1)
1008 interiorElementBoundaries.insert(ebN);
1012 exteriorElementBoundaries.insert(ebN);
1014 if (eb->second.right != -1)
1024 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1026 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1028 set<NodeTuple<2> > edges;
1043 for (set<
NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
1060 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1062 stop = CurrentTime();
1080 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1123 using namespace std;
1127 ElementBoundaryElementsType elementBoundaryElements;
1128 start=CurrentTime();
1138 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1140 elementBoundaryElements[ebt].right=eN;
1141 elementBoundaryElements[ebt].right_ebN_element=ebN;
1145 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1148 stop = CurrentTime();
1154 start = CurrentTime();
1155 typedef set<int> ElementBoundariesType;
1156 ElementBoundariesType interiorElementBoundaries,exteriorElementBoundaries;
1162 stop = CurrentTime();
1166 start = CurrentTime();
1169 eb != elementBoundaryElements.end();
1181 if(eb->second.right != -1)
1183 interiorElementBoundaries.insert(ebN);
1187 exteriorElementBoundaries.insert(ebN);
1189 if (eb->second.right != -1)
1199 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1201 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1208 typedef set<NodeTuple<2> > EdgesType;
1224 EdgesType::iterator edge_p=edges.begin();
1225 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
1247 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1252 stop = CurrentTime();
1270 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1312 using namespace std;
1316 start=CurrentTime();
1318 int lface[6][4] = {{0,1,2,3},
1336 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1338 elementBoundaryElements[ebt].right=eN;
1339 elementBoundaryElements[ebt].right_ebN_element=ebN;
1343 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1346 stop = CurrentTime();
1352 start = CurrentTime();
1353 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1360 stop = CurrentTime();
1364 start = CurrentTime();
1367 eb != elementBoundaryElements.end();
1380 if(eb->second.right != -1)
1382 interiorElementBoundaries.insert(ebN);
1386 exteriorElementBoundaries.insert(ebN);
1389 if (eb->second.right != -1)
1399 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1401 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1403 set<NodeTuple<2> > edges;
1405 int ledge[12][2] = {{0,1},{1,2},{2,3},{3,0},
1406 {0,4},{1,5},{2,6},{3,7},
1407 {4,5},{5,6},{6,7},{7,4}};
1416 for (
int e=0;e<12;e++)
1429 set<NodeTuple<2> >::iterator edge_p=edges.begin();
1430 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
1453 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1455 stop = CurrentTime();
1475 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1518 using namespace std;
1522 int n2 = (mesh.
px+1)*(mesh.
py+1)-1 ;
1523 int n3 = (mesh.
px+1)*mesh.
py;
1526 int n4 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n0;
1527 int n5 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n1;
1528 int n6 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n2;
1529 int n7 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n3;
1535 using namespace std;
1539 start=CurrentTime();
1541 int lface[6][4] = {{n0,n1,n2,n3},
1561 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1563 elementBoundaryElements[ebt].right=eN;
1564 elementBoundaryElements[ebt].right_ebN_element=ebN;
1568 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1571 stop = CurrentTime();
1577 start = CurrentTime();
1578 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1585 stop = CurrentTime();
1589 start = CurrentTime();
1592 eb != elementBoundaryElements.end();
1605 if(eb->second.right != -1)
1607 interiorElementBoundaries.insert(ebN);
1611 exteriorElementBoundaries.insert(ebN);
1614 if (eb->second.right != -1)
1624 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1626 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1628 set<NodeTuple<2> > edges;
1633 int ledge[12][2] = {{n0,n1},{n1,n2},{n2,n3},{n3,n0},
1634 {n0,n4},{n1,n5},{n2,n6},{n3,n7},
1635 {n4,n5},{n5,n6},{n6,n7},{n7,n4}};
1644 for (
int e=0;e<12;e++)
1657 set<NodeTuple<2> >::iterator edge_p=edges.begin();
1658 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
1678 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1680 stop = CurrentTime();
1699 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1743 int n2 = (mesh.
px+1)*(mesh.
py+1)-1 ;
1744 int n3 = (mesh.
px+1)*mesh.
py;
1747 int n4 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n0;
1748 int n5 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n1;
1749 int n6 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n2;
1750 int n7 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n3;
1755 using namespace std;
1760 int> elementBoundaryIds;
1761 start=CurrentTime();
1763 int lface[6][4] = {{n0,n1,n2,n3},
1781 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1783 elementBoundaryElements[ebt].right=eN;
1784 elementBoundaryElements[ebt].right_ebN_element=ebN;
1790 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1791 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
1794 stop = CurrentTime();
1800 start = CurrentTime();
1801 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1806 stop = CurrentTime();
1810 start = CurrentTime();
1812 eb != elementBoundaryElements.end();
1815 int ebN = elementBoundaryIds[eb->first];
1826 if(eb->second.right != -1)
1828 interiorElementBoundaries.insert(ebN);
1832 exteriorElementBoundaries.insert(ebN);
1834 if (eb->second.right != -1)
1844 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1846 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1886 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1888 stop = CurrentTime();
1906 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1949 using namespace std;
1955 int> elementBoundaryIds;
1956 start=CurrentTime();
1958 for(
int ebN=0;ebN<2;ebN++)
1964 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1966 elementBoundaryElements[ebt].right=eN;
1967 elementBoundaryElements[ebt].right_ebN_element=ebN;
1968 assert(elementBoundaryIds[ebt] == ebN_global);
1972 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1973 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
1976 stop = CurrentTime();
1982 start = CurrentTime();
1983 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1988 stop = CurrentTime();
1992 start = CurrentTime();
1994 eb != elementBoundaryElements.end();
1997 int ebN = elementBoundaryIds[eb->first];
2004 if(eb->second.right != -1)
2006 interiorElementBoundaries.insert(ebN);
2010 exteriorElementBoundaries.insert(ebN);
2012 if (eb->second.right != -1)
2022 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2024 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2046 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2051 stop = CurrentTime();
2066 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2107 using namespace std;
2116 int> elementBoundaryIds;
2117 start=CurrentTime();
2119 for(
int ebN=0;ebN<3;ebN++)
2126 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2128 elementBoundaryElements[ebt].right=eN;
2129 elementBoundaryElements[ebt].right_ebN_element=ebN;
2130 assert(elementBoundaryIds[ebt] == ebN_global);
2134 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2135 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2138 stop = CurrentTime();
2144 start = CurrentTime();
2145 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2150 stop = CurrentTime();
2154 start = CurrentTime();
2156 eb != elementBoundaryElements.end();
2159 int ebN = elementBoundaryIds[eb->first];
2167 if(eb->second.right != -1)
2169 interiorElementBoundaries.insert(ebN);
2173 exteriorElementBoundaries.insert(ebN);
2175 if (eb->second.right != -1)
2185 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2187 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2189 set<NodeTuple<2> > edges;
2204 for (set<
NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
2221 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2223 stop = CurrentTime();
2241 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2281 using namespace std;
2290 int> elementBoundaryIds;
2291 start=CurrentTime();
2293 for(
int ebN=0;ebN<4;ebN++)
2300 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2302 elementBoundaryElements[ebt].right=eN;
2303 elementBoundaryElements[ebt].right_ebN_element=ebN;
2304 assert(elementBoundaryIds[ebt] == ebN_global);
2308 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2309 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2312 stop = CurrentTime();
2318 start = CurrentTime();
2319 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2324 stop = CurrentTime();
2328 start = CurrentTime();
2330 eb != elementBoundaryElements.end();
2333 int ebN = elementBoundaryIds[eb->first];
2341 if(eb->second.right != -1)
2343 interiorElementBoundaries.insert(ebN);
2347 exteriorElementBoundaries.insert(ebN);
2349 if (eb->second.right != -1)
2359 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2361 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2363 set<NodeTuple<2> > edges;
2377 for (set<
NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
2394 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2396 stop = CurrentTime();
2414 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2457 using namespace std;
2462 int> elementBoundaryIds;
2463 start=CurrentTime();
2474 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2476 elementBoundaryElements[ebt].right=eN;
2477 elementBoundaryElements[ebt].right_ebN_element=ebN;
2478 assert(elementBoundaryIds[ebt] == ebN_global);
2482 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2483 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2486 stop = CurrentTime();
2492 start = CurrentTime();
2493 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2498 stop = CurrentTime();
2502 start = CurrentTime();
2504 eb != elementBoundaryElements.end();
2507 int ebN = elementBoundaryIds[eb->first];
2517 if(eb->second.right != -1)
2519 interiorElementBoundaries.insert(ebN);
2523 exteriorElementBoundaries.insert(ebN);
2525 if (eb->second.right != -1)
2535 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2537 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2539 set<NodeTuple<2> > edges;
2554 set<NodeTuple<2> >::iterator edge_p=edges.begin();
2555 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
2573 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2575 stop = CurrentTime();
2593 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2637 using namespace std;
2643 int> elementBoundaryIds;
2644 start=CurrentTime();
2646 for(
int ebN=0;ebN<2;ebN++)
2652 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2654 elementBoundaryElements[ebt].right=eN;
2655 elementBoundaryElements[ebt].right_ebN_element=ebN;
2656 assert(elementBoundaryIds[ebt] == ebN_global);
2660 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2661 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2664 stop = CurrentTime();
2670 start = CurrentTime();
2671 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2676 stop = CurrentTime();
2680 start = CurrentTime();
2682 eb != elementBoundaryElements.end();
2685 int ebN = elementBoundaryIds[eb->first];
2692 if(eb->second.right != -1)
2694 interiorElementBoundaries.insert(ebN);
2698 exteriorElementBoundaries.insert(ebN);
2700 if (eb->second.right != -1)
2710 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2712 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2728 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2733 stop = CurrentTime();
2748 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2789 using namespace std;
2798 int> elementBoundaryIds;
2799 start=CurrentTime();
2801 for(
int ebN=0;ebN<3;ebN++)
2808 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2810 elementBoundaryElements[ebt].right=eN;
2811 elementBoundaryElements[ebt].right_ebN_element=ebN;
2812 assert(elementBoundaryIds[ebt] == ebN_global);
2816 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2817 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2820 stop = CurrentTime();
2826 start = CurrentTime();
2827 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2832 stop = CurrentTime();
2836 start = CurrentTime();
2838 eb != elementBoundaryElements.end();
2841 int ebN = elementBoundaryIds[eb->first];
2849 if(eb->second.right != -1)
2851 interiorElementBoundaries.insert(ebN);
2855 exteriorElementBoundaries.insert(ebN);
2857 if (eb->second.right != -1)
2867 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2869 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2904 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2906 stop = CurrentTime();
2924 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2964 using namespace std;
2973 int> elementBoundaryIds;
2974 start=CurrentTime();
2976 for(
int ebN=0;ebN<4;ebN++)
2983 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2985 elementBoundaryElements[ebt].right=eN;
2986 elementBoundaryElements[ebt].right_ebN_element=ebN;
2987 assert(elementBoundaryIds[ebt] == ebN_global);
2991 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2992 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2995 stop = CurrentTime();
3001 start = CurrentTime();
3002 set<int> interiorElementBoundaries,exteriorElementBoundaries;
3007 stop = CurrentTime();
3011 start = CurrentTime();
3013 eb != elementBoundaryElements.end();
3016 int ebN = elementBoundaryIds[eb->first];
3024 if(eb->second.right != -1)
3026 interiorElementBoundaries.insert(ebN);
3030 exteriorElementBoundaries.insert(ebN);
3032 if (eb->second.right != -1)
3042 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
3044 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
3059 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3061 stop = CurrentTime();
3079 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3129 using namespace std;
3134 int> elementBoundaryIds;
3135 start=CurrentTime();
3148 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
3150 elementBoundaryElements[ebt].right=eN;
3151 elementBoundaryElements[ebt].right_ebN_element=ebN;
3152 assert(elementBoundaryIds[ebt] == ebN_global);
3156 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
3157 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
3160 stop = CurrentTime();
3166 start = CurrentTime();
3167 set<int> interiorElementBoundaries,exteriorElementBoundaries;
3172 stop = CurrentTime();
3176 start = CurrentTime();
3178 eb != elementBoundaryElements.end();
3182 int ebN = elementBoundaryIds[eb->first];
3192 if(eb->second.right != -1)
3194 interiorElementBoundaries.insert(ebN);
3198 exteriorElementBoundaries.insert(ebN);
3200 if (eb->second.right != -1)
3210 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
3212 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
3249 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3268 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3310 using namespace std;
3315 int> elementBoundaryIds;
3316 start=CurrentTime();
3318 int lface[6][4] = {{0,1,2,3},
3336 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
3338 elementBoundaryElements[ebt].right=eN;
3339 elementBoundaryElements[ebt].right_ebN_element=ebN;
3343 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
3344 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
3347 stop = CurrentTime();
3353 start = CurrentTime();
3354 set<int> interiorElementBoundaries,exteriorElementBoundaries;
3359 stop = CurrentTime();
3363 start = CurrentTime();
3365 eb != elementBoundaryElements.end();
3368 int ebN = elementBoundaryIds[eb->first];
3379 if(eb->second.right != -1)
3381 interiorElementBoundaries.insert(ebN);
3385 exteriorElementBoundaries.insert(ebN);
3387 if (eb->second.right != -1)
3397 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
3399 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
3437 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3439 stop = CurrentTime();
3457 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3498 dx = nodeArray[nL*3+0] - nodeArray[nR*3+0];
3499 dy = nodeArray[nL*3+1] - nodeArray[nR*3+1];
3500 dz = nodeArray[nL*3+2] - nodeArray[nR*3+2];
3501 return sqrt(dx*dx+dy*dy+dz*dz);
3504 inline double triangleArea(
int n0,
int n1,
int n2,
const double* nodeArray)
3507 va[0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
3508 va[1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
3509 va[2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
3510 vb[0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
3511 vb[1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
3512 vb[2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
3513 return 0.5*fabs(va[1]*vb[2] - vb[1]*va[2] - (va[0]*vb[2] - vb[0]*va[2]) + (va[0]*vb[1] - va[1]*vb[0]));
3519 t[0][0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
3520 t[0][1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
3521 t[0][2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
3523 t[1][0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
3524 t[1][1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
3525 t[1][2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
3527 t[2][0] = nodeArray[n3*3+0] - nodeArray[n0*3+0];
3528 t[2][1] = nodeArray[n3*3+1] - nodeArray[n0*3+1];
3529 t[2][2] = nodeArray[n3*3+2] - nodeArray[n0*3+2];
3530 return fabs(t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) - \
3531 t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) + \
3532 t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]))/6.0;
3588 double volume, surfaceArea=0.0, hMin=0.0,hMax=0.0;
3595 for (
int ebN=0;ebN<4;ebN++)
3602 hMin = 6.0*volume/surfaceArea;
3614 mesh.
h = fmax(hMax,mesh.
h);
3642 inline double hexahedronVolume(
int n0,
int n1,
int n2,
int n3,
int n4,
int n5,
int n6,
int n7,
const double* nodeArray)
3645 t[0] = nodeArray[n0*3+0] - nodeArray[n1*3+0];
3646 t[1] = nodeArray[n0*3+1] - nodeArray[n3*3+1];
3647 t[2] = nodeArray[n0*3+2] - nodeArray[n4*3+2];
3649 return fabs(t[0]*t[1]*t[2]);
3692 double volume, hMin=0.0,hMax=0.0;
3722 mesh.
h = fmax(hMax,mesh.
h);
3789 double area=0.0, perimeter=0.0,h=0.0, hMin=0.0,hMax=0.0;
3802 hMax = fmax(hMax,h);
3804 hMin = 4.0*area/perimeter;
3809 mesh.
h = fmax(hMax,mesh.
h);
3885 double area=0.0, perimeter=0.0,h=0.0, hMin=1.0e6,hMax=0.0;
3903 hMax = fmax(hMax,h);
3904 hMin = fmin(hMin,h);
3910 mesh.
h = fmax(hMax,mesh.
h);
4029 using namespace std;
4030 multilevelMesh.
nLevels = nLevels;
4032 for(
int i=0;i<nLevels;i++)
4038 for(
int i=1;i<nLevels;i++)
4041 set<Node> newNodeSet;
4042 set<Node>::iterator nodeItr;
4043 pair<set<Node>::iterator,
bool> ret;
4067 midpoints[0].
nN = nN_new;
4071 else if (averageNewNodeFlags)
4077 newNodeSet.insert(midpoints[0]);
4097 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4114 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4124 using namespace std;
4125 multilevelMesh.
nLevels = nLevels;
4127 for(
int i=0;i<nLevels;i++)
4133 for(
int i=1;i<nLevels;i++)
4136 set<Node> newNodeSet;
4137 set<Node>::iterator nodeItr;
4138 pair<set<Node>::iterator,
bool> ret;
4165 for(
int nN_element_0=0,nN_midpoint=0;nN_element_0<mesh.
nNodes_element;nN_element_0++)
4166 for(
int nN_element_1=nN_element_0+1;nN_element_1<mesh.
nNodes_element;nN_element_1++,nN_midpoint++)
4170 midpoints[nN_midpoint]);
4171 nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
4172 if(nodeItr == newNodeSet.end())
4174 midpoints[nN_midpoint].
nN = nN_new;
4178 else if (averageNewNodeFlags)
4184 newNodeSet.insert(midpoints[nN_midpoint]);
4188 midpoints[nN_midpoint].
nN = nodeItr->nN;
4218 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4235 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4249 using namespace std;
4250 multilevelMesh.
nLevels = nLevels;
4252 for(
int i=0;i<nLevels;i++)
4258 for(
int i=1;i<nLevels;i++)
4260 std::cout<<
"Quad refinement not imlemented "<<i<<endl;
4268 using namespace std;
4269 multilevelMesh.
nLevels = nLevels;
4271 for(
int i=0;i<nLevels;i++)
4277 for(
int i=1;i<nLevels;i++)
4280 set<Node> newNodeSet;
4281 set<Node>::iterator nodeItr;
4282 pair<set<Node>::iterator,
bool> ret;
4324 for(
int nN_element_0=0,nN_midpoint=0;nN_element_0<mesh.
nNodes_element;nN_element_0++)
4325 for(
int nN_element_1=nN_element_0+1;nN_element_1<mesh.
nNodes_element;nN_element_1++,nN_midpoint++)
4329 midpoints[nN_midpoint]);
4330 nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
4331 if(nodeItr == newNodeSet.end())
4333 midpoints[nN_midpoint].
nN = nN_new;
4337 else if (averageNewNodeFlags)
4343 newNodeSet.insert(midpoints[nN_midpoint]);
4347 midpoints[nN_midpoint].
nN = nodeItr->nN;
4373 for(
int dN=1;dN<3;dN++)
4375 dd =
edgeLength(midpoints[dN],midpoints[5-dN]);
4381 else if (dd == mind)
4383 if(midpoints[dN] < midpoints[mindN])
4385 else if (!(midpoints[mindN] < midpoints[dN]) && (midpoints[5-dN] < midpoints[5-mindN]))
4412 else if (mindN == 1)
4469 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4543 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4553 using namespace std;
4554 multilevelMesh.
nLevels = nLevels;
4556 for(
int i=0;i<nLevels;i++)
4562 for(
int i=1;i<nLevels;i++)
4564 std::cout<<
"Hexahedron refinement not implemented"<<std::endl;
4571 const int& nSpace_global)
4579 int eN_left_parent = levelElementParentsArray[eN_left];
4580 int eN_right_parent= levelElementParentsArray[eN_right];
4581 if (eN_left_parent == eN_right_parent)
4586 int left_parent_ebN = -1;
4588 ebN_left_parent_element++)
4591 ebN_left_parent_element]
4594 left_parent_ebN = ebN_left_parent_element;
4608 int eN_parent = levelElementParentsArray[eN];
4609 int lastExteriorElementBoundaryOnParent=-1;
4610 int nExteriorElementBoundariesOnParent =0;
4615 lastExteriorElementBoundaryOnParent = ebN_element;
4616 nExteriorElementBoundariesOnParent++;
4619 assert(nExteriorElementBoundariesOnParent > 0); assert(lastExteriorElementBoundaryOnParent >= 0);
4620 if (nExteriorElementBoundariesOnParent > 1)
4624 if (nSpace_global == 1)
4631 const double n_child = (childMesh.
nodeArray[nN_child1*3+0]-childMesh.
nodeArray[nN_child0*3+0])/
4636 const double n_parent0 = (parentMesh.
nodeArray[nN_parent1*3+0]-parentMesh.
nodeArray[nN_parent0*3+0])/
4638 const double n_parent1 = (parentMesh.
nodeArray[nN_parent0*3+0]-parentMesh.
nodeArray[nN_parent1*3+0])/
4641 int ebN_parent = -1;
4642 if (fabs(n_child-n_parent1) < 1.0e-8)
4644 else if (fabs(n_child-n_parent0) < 1.0e-8)
4646 assert(ebN_parent >= 0);
4649 else if (nSpace_global == 2)
4653 double n_child[2] = {childMesh.
nodeArray[nN_ebN_child[1]*3+1]-childMesh.
nodeArray[nN_ebN_child[0]*3+1],
4655 double tmp = sqrt(n_child[0]*n_child[0] + n_child[1]*n_child[1]);
4657 n_child[0] /= tmp; n_child[1] /= tmp;
4659 int ebN_parent = -1;
4666 double n_parent[2] = {parentMesh.
nodeArray[nN_ebN_parent[1]*3+1]-parentMesh.
nodeArray[nN_ebN_parent[0]*3+1],
4667 parentMesh.
nodeArray[nN_ebN_parent[0]*3+0]-parentMesh.
nodeArray[nN_ebN_parent[1]*3+0]};
4668 tmp = sqrt(n_parent[0]*n_parent[0] + n_parent[1]*n_parent[1]);
4670 n_parent[0] /= tmp; n_parent[1] /= tmp;
4673 tmp = n_parent[0]*n_child[0] + n_parent[1]*n_child[1];
4675 if (fabs(sqrt(tmp*tmp) - 1.0) < 1.0e-8)
4677 ebN_parent = ebN_cur;
4681 assert(ebN_parent >= 0);
4689 double n_child[3] = {0.,0.,0.};
4691 n_child[0] = (nN_ebN_child_1_x[1]-nN_ebN_child_0_x[1])*(nN_ebN_child_2_x[2]-nN_ebN_child_0_x[2])
4692 -(nN_ebN_child_1_x[2]-nN_ebN_child_0_x[2])*(nN_ebN_child_2_x[1]-nN_ebN_child_0_x[1]);
4694 n_child[1] = (nN_ebN_child_1_x[2]-nN_ebN_child_0_x[2])*(nN_ebN_child_2_x[0]-nN_ebN_child_0_x[0])
4695 -(nN_ebN_child_1_x[0]-nN_ebN_child_0_x[0])*(nN_ebN_child_2_x[2]-nN_ebN_child_0_x[2]);
4697 n_child[2] = (nN_ebN_child_1_x[0]-nN_ebN_child_0_x[0])*(nN_ebN_child_2_x[1]-nN_ebN_child_0_x[1])
4698 -(nN_ebN_child_1_x[1]-nN_ebN_child_0_x[1])*(nN_ebN_child_2_x[0]-nN_ebN_child_0_x[0]);
4700 double tmp = sqrt(n_child[0]*n_child[0] + n_child[1]*n_child[1] + n_child[2]*n_child[2]);
4702 n_child[0] /= tmp; n_child[1] /= tmp; n_child[2] /= tmp;
4704 int ebN_parent = -1;
4705 double n_parent[3] = {0.,0.,0.};
4714 n_parent[0] = (nN_ebN_parent_1_x[1]-nN_ebN_parent_0_x[1])*(nN_ebN_parent_2_x[2]-nN_ebN_parent_0_x[2])
4715 -(nN_ebN_parent_1_x[2]-nN_ebN_parent_0_x[2])*(nN_ebN_parent_2_x[1]-nN_ebN_parent_0_x[1]);
4717 n_parent[1] = (nN_ebN_parent_1_x[2]-nN_ebN_parent_0_x[2])*(nN_ebN_parent_2_x[0]-nN_ebN_parent_0_x[0])
4718 -(nN_ebN_parent_1_x[0]-nN_ebN_parent_0_x[0])*(nN_ebN_parent_2_x[2]-nN_ebN_parent_0_x[2]);
4720 n_parent[2] = (nN_ebN_parent_1_x[0]-nN_ebN_parent_0_x[0])*(nN_ebN_parent_2_x[1]-nN_ebN_parent_0_x[1])
4721 -(nN_ebN_parent_1_x[1]-nN_ebN_parent_0_x[1])*(nN_ebN_parent_2_x[0]-nN_ebN_parent_0_x[0]);
4723 tmp = sqrt(n_parent[0]*n_parent[0] + n_parent[1]*n_parent[1] + n_parent[2]*n_parent[2]);
4725 n_parent[0] /= tmp; n_parent[1] /= tmp; n_parent[2] /= tmp;
4729 tmp = n_parent[0]*n_child[0] + n_parent[1]*n_child[1] + n_parent[2]*n_child[2];
4731 if (fabs(sqrt(tmp*tmp) - 1.0) < 1.0e-8)
4733 ebN_parent = ebN_cur;
4738 assert(ebN_parent >= 0);
4746 lastExteriorElementBoundaryOnParent];
4758 using namespace std;
4761 string word,elementType;
4766 if(word ==
"MESH1D")
4768 elementType =
"E2E";
4772 else if(word ==
"MESH2D")
4774 elementType =
"E3T";
4778 else if (word ==
"MESH3D")
4780 elementType =
"E4T";
4786 cerr<<
"Unrecognized mesh type"<<endl;
4794 while(word == elementType)
4799 meshFile>>nodes[nN];
4802 elementNodesVector.push_back(nodes);
4804 materialTypes.push_back(type);
4822 using namespace std;
4824 meshFile<<
"ND"<<setw(7)<<nN+1
4825 <<scientific<<setprecision(8)<<setw(16)<<mesh.
nodeArray[3*nN+0]
4826 <<scientific<<setprecision(8)<<setw(16)<<mesh.
nodeArray[3*nN+1]
4827 <<scientific<<setprecision(8)<<setw(16)<<mesh.
nodeArray[3*nN+2]
4835 meshFile<<
"Try to write something"<<std::endl;
4836 using namespace std;
4842 elementType =
"E2E";
4843 meshFile<<
"MESH1D"<<endl;
4848 elementType =
"E3T";
4849 meshFile<<
"MESH2D"<<endl;
4854 elementType =
"E4T";
4855 meshFile<<
"MESH3D"<<endl;
4859 cerr<<
"Unknown element type"<<endl;
4864 meshFile<<elementType;
4865 meshFile<<setw(width)<<eN+1;
4883 using namespace std;
4887 int nElements_file, nNodes_file;
4897 elementNodesArray_file,
4898 nodeMaterialTypes_file,
4899 elementMaterialTypes_file,
4915 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
4918 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
4927 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
4929 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
4941 using namespace std;
4975 using namespace std;
4979 int nElementBoundaries_file;
4982 bool elementBoundaryMaterialTypesInFile =
false;
4985 elementBoundaryMaterialTypesInFile,
4986 nElementBoundaries_file,
4987 elementBoundaryNodesArray_file,
4988 elementBoundaryMaterialTypes_file,
4995 if (!elementBoundaryMaterialTypesInFile)
5001 map<NodeTuple<2>,
int> triangleElementBoundaryMaterialTypes;
5003 for (
int ebN = 0; ebN < nElementBoundaries_file; ebN++)
5006 nodes[0] = elementBoundaryNodesArray_file[ebN*2+0];
5007 nodes[1] = elementBoundaryNodesArray_file[ebN*2+1];
5009 triangleElementBoundaryMaterialTypes[ttuple] = elementBoundaryMaterialTypes_file[ebN];
5031 using namespace std;
5035 int nElements_file, nNodes_file;
5045 elementNodesArray_file,
5046 nodeMaterialTypes_file,
5047 elementMaterialTypes_file,
5062 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5066 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5075 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5077 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5092 using namespace std;
5096 int nElementBoundaries_file;
5099 bool elementBoundaryMaterialTypesInFile =
false;
5102 elementBoundaryMaterialTypesInFile,
5103 nElementBoundaries_file,
5104 elementBoundaryNodesArray_file,
5105 elementBoundaryMaterialTypes_file,
5112 if (!elementBoundaryMaterialTypesInFile)
5118 using namespace std;
5123 start=CurrentTime();
5133 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
5135 elementBoundaryElements[ebt].right=eN;
5136 elementBoundaryElements[ebt].right_ebN_element=ebN;
5140 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
5143 stop = CurrentTime();
5149 start = CurrentTime();
5157 stop = CurrentTime();
5161 start = CurrentTime();
5164 eb != elementBoundaryElements.end();
5185 if (eb->second.right != -1)
5202 set<NodeTuple<2> > edges;
5217 set<NodeTuple<2> >::iterator edge_p=edges.begin();
5218 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
5237 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
5239 stop = CurrentTime();
5258 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
5306 map<NodeTuple<3>,
int> tetgenElementBoundaryMaterialTypes;
5308 for (
int ebN = 0; ebN < nElementBoundaries_file; ebN++)
5311 nodes[0] = elementBoundaryNodesArray_file[ebN*3+0];
5312 nodes[1] = elementBoundaryNodesArray_file[ebN*3+1];
5313 nodes[2] = elementBoundaryNodesArray_file[ebN*3+2];
5315 tetgenElementBoundaryMaterialTypes[ttuple] = elementBoundaryMaterialTypes_file[ebN];
5332 map<NodeTuple<3>,
int> tetgenElementBoundaryMaterialTypes;
5334 for (
int ebNE = 0; ebNE < nElementBoundaries_file; ebNE++)
5337 nodes[0] = elementBoundaryNodesArray_file[ebNE*3+0];
5338 nodes[1] = elementBoundaryNodesArray_file[ebNE*3+1];
5339 nodes[2] = elementBoundaryNodesArray_file[ebNE*3+2];
5341 tetgenElementBoundaryMaterialTypes[ttuple] = elementBoundaryMaterialTypes_file[ebNE];
5356 printf(
"Assertion failed in readTetgenElementBoundaryMaterialTypes\n");
5376 using namespace std;
5393 bool writeExteriorElementBoundariesOnly =
true;
5397 nElementBoundariesToWrite,
5400 writeExteriorElementBoundariesOnly,
5413 using namespace std;
5416 std::string meshFilename= std::string(filebase)+
".3dm";
5417 std::ifstream meshFile(meshFilename.c_str());
5418 if (!meshFile.good())
5420 std::cerr<<
"read3DM cannot open file "
5421 <<meshFilename<<std::endl;
5426 std::string fileType;
5428 if (fileType !=
"MESH3D")
5430 std::cerr<<
"read3DM does not recognize filetype "
5431 <<fileType<<std::endl;
5435 std::string firstWord;
5436 meshFile>>firstWord;
5438 std::vector<int> elementNodesArray_file;
5439 std::vector<int> elementMaterialTypes_file;
5440 int eN,n0,n1,n2,n3,emt;
5441 while (firstWord ==
"E4T")
5444 meshFile>>eN>>n0>>n1>>n2>>n3>>emt;
5445 elementNodesArray_file.push_back(n0-indexBase);
5446 elementNodesArray_file.push_back(n1-indexBase);
5447 elementNodesArray_file.push_back(n2-indexBase);
5448 elementNodesArray_file.push_back(n3-indexBase);
5449 elementMaterialTypes_file.push_back(emt-indexBase);
5450 meshFile>>firstWord;
5452 std::vector<int> nodeMaterialTypes_file;
5453 std::vector<double> nodeArray_file;
5457 while (!meshFile.eof() && firstWord ==
"ND")
5460 meshFile>>nN>>x>>y>>
z;
5461 nodeArray_file.push_back(x);
5462 nodeArray_file.push_back(y);
5463 nodeArray_file.push_back(
z);
5464 nodeMaterialTypes_file.push_back(0);
5465 meshFile>>firstWord;
5470 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5474 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5479 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5484 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5497 using namespace std;
5500 std::string meshFilename= std::string(filebase)+
".3dm";
5501 std::ifstream meshFile(meshFilename.c_str());
5502 if (!meshFile.good())
5504 std::cerr<<
"read2DM cannot open file "
5505 <<meshFilename<<std::endl;
5510 std::string fileType;
5512 if (fileType !=
"MESH2D")
5514 std::cerr<<
"read2DM does not recognize filetype "
5515 <<fileType<<std::endl;
5519 std::string meshName;
5520 std::string firstWord;
5522 if (meshName ==
"E3T")
5523 firstWord = meshName;
5527 logEvent(&(
"Reading 2DM"+meshName)[0],5);
5528 meshFile>>firstWord;
5531 std::vector<int> elementNodesArray_file;
5532 std::vector<int> elementMaterialTypes_file;
5533 int eN,n0,n1,n2,emt;
5534 while (firstWord ==
"E3T")
5537 meshFile>>eN>>n0>>n1>>n2>>emt;
5538 elementNodesArray_file.push_back(n0-indexBase);
5539 elementNodesArray_file.push_back(n1-indexBase);
5540 elementNodesArray_file.push_back(n2-indexBase);
5541 elementMaterialTypes_file.push_back(emt-indexBase);
5542 meshFile>>firstWord;
5544 std::vector<int> nodeMaterialTypes_file;
5545 std::vector<double> nodeArray_file;
5549 while (!meshFile.eof() && firstWord ==
"ND")
5552 meshFile>>nN>>x>>y>>
z;
5553 nodeArray_file.push_back(x);
5554 nodeArray_file.push_back(y);
5555 nodeArray_file.push_back(
z);
5556 nodeMaterialTypes_file.push_back(0);
5557 meshFile>>firstWord;
5562 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5566 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5571 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5576 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5589 using namespace std;
5592 std::string meshFilename= std::string(filebase)+
".mesh";
5593 std::ifstream meshFile(meshFilename.c_str());
5597 if (!meshFile.good())
5599 std::cerr<<
"readHex cannot open file "
5600 <<meshFilename<<std::endl;
5605 std::string fileType;
5607 if (fileType !=
"HEX")
5609 std::cerr<<
"readHex does not recognize filetype "
5610 <<fileType<<std::endl;
5635 int n0,n1,n2,n3,n4,n5,n6,n7,emt;
5639 meshFile>>n0>>n1>>n2>>n3>>n4>>n5>>n6>>n7>>emt;
5664 using namespace std;
5667 std::string bcFilename= std::string(filebase)+
".bc";
5668 std::ifstream bcFile(bcFilename.c_str());
5671 std::cerr<<
"readBC cannot open file "
5672 <<bcFilename<<std::endl;
5676 std::string firstWord;
5677 int eN,ebN_local,nN,flag;
5678 while (!bcFile.eof())
5681 if (firstWord ==
"FCS")
5683 bcFile>>eN>>ebN_local>>flag;
5686 if (firstWord ==
"NDS")
5700 using namespace std;
5720 using namespace std;
5739 using namespace std;
5742 int nLevelsNew = multilevelMesh.
nLevels+nLevels2add;
5743 Mesh * meshArrayTmp =
new Mesh[nLevelsNew];
5744 int** elementChildrenArrayTmp =
new int*[nLevelsNew];
5745 int** elementChildrenOffsetsTmp =
new int*[nLevelsNew];
5746 int** elementParentsArrayTmp =
new int*[nLevelsNew];
5749 for (
int i=0; i < multilevelMesh.
nLevels; i++)
5751 meshArrayTmp[i] = multilevelMesh.
meshArray[i];
5761 multilevelMesh.
meshArray = meshArrayTmp;
5767 for (
int i=multilevelMesh.
nLevels; i < nLevelsNew; i++)
5774 multilevelMesh.
nLevels = nLevelsNew;
5779 int * elementTagArray)
5781 using namespace std;
5787 int nLevelsPrev = multilevelMesh.
nLevels;
5789 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
5794 int nElements_tagged = 0;
5797 if (elementTagArray[eN] > 0)
5851 set<Node> newNodeSet;
5852 set<Node>::iterator nodeItr;
5857 if (elementTagArray[eN_parent] == 0)
5871 else if (elementTagArray[eN_parent] > 0)
5888 midpoints[0].
nN = nN_new;
5889 newNodeSet.insert(midpoints[0]);
5906 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
5929 int * elementTagArray)
5931 using namespace std;
5937 int nLevelsPrev = multilevelMesh.
nLevels;
5941 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
5942 int nElements_tagged = 0;
5945 if (elementTagArray[eN] > 0)
5951 vector<int> elementNodesArray_tmp,elementNeighborsArray_tmp,bases_tmp,
5952 elementParentsArray_tmp;
5960 4*nElements_tagged));
5962 4*nElements_tagged));
5964 4*nElements_tagged);
5967 3*nElements_tagged));
5971 4*nElements_tagged);
5975 elementNodesArray_tmp.insert(elementNodesArray_tmp.begin(),
5980 elementNeighborsArray_tmp.insert(elementNeighborsArray_tmp.begin(),
5985 bases_tmp.insert(bases_tmp.begin(),
5991 nodeArray_tmp.insert(nodeArray_tmp.begin(),
6002 elementParentsArray_tmp.push_back(eN);
6004 for (
int eN_parent = 0;
6009 if (elementTagArray[eN_parent] == 1 && !refined[eN_parent])
6015 elementNodesArray_tmp,
6016 elementNeighborsArray_tmp,
6017 elementChildrenList,
6018 elementParentsArray_tmp,
6025 assert(elementParentsArray_tmp.size() ==
unsigned(nElements_new));
6028 copy(elementParentsArray_tmp.begin(),elementParentsArray_tmp.end(),
6039 if (elementChildrenList[eN].size() > 0)
6042 offset + elementChildrenList[eN].size();
6044 for (std::list<int>::iterator it = elementChildrenList[eN].begin();
6045 it != elementChildrenList[eN].end(); it++)
6080 copy(elementNodesArray_tmp.begin(),elementNodesArray_tmp.end(),
6082 copy(nodeArray_tmp.begin(),nodeArray_tmp.end(),
6084 copy(bases_tmp.begin(),bases_tmp.end(),
6109 using namespace std;
6110 int nLevels = multilevelMesh.
nLevels;
6113 logEvent(
"WARNING setNewestNodeBasesToLongestEdge newestNodeBases !=NULL exiting",5);
6119 logEvent(
"WARNING setNewestNodeBasesToLongestEdge 2d only exiting",5);
6136 int * elementTagArray)
6138 using namespace std;
6141 int nLevelsPrev = multilevelMesh.
nLevels;
6143 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
6146 set<Node> newNodeSet;
6147 set<Node>::iterator nodeItr;
6149 map<int, vector<int> > elementChildren;
6150 map<int, int> elementParents;
6152 map<int,int> elementsForBisection;
6153 map<int,vector<int> > elementsForTrisection,newElementsForTrisection,nextElementsForTrisection;
6154 set<int> elementsForUniform,newElementsForUniform,nextElementsForUniform;
6158 set<int> oldElements;
6159 int i = nLevelsPrev;
6163 oldElements.insert(eN_parent);
6164 if (elementTagArray[eN_parent] > 0)
6166 newElementsForUniform.insert(eN_parent);
6167 elementsForUniform.insert(eN_parent);
6172 while (!newElementsForUniform.empty() || !newElementsForTrisection.empty())
6175 for(set<int>::iterator eN_uniform_itr = newElementsForUniform.begin(); eN_uniform_itr != newElementsForUniform.end(); eN_uniform_itr++)
6177 int eN_parent = *eN_uniform_itr;
6178 oldElements.erase(eN_parent);
6179 elementChildren[eN_parent].push_back(eN_new+0);
6180 elementChildren[eN_parent].push_back(eN_new+1);
6181 elementChildren[eN_parent].push_back(eN_new+2);
6182 elementChildren[eN_parent].push_back(eN_parent);
6183 elementParents[eN_new + 0] = eN_parent;
6184 elementParents[eN_new + 1] = eN_parent;
6185 elementParents[eN_new + 2] = eN_parent;
6186 elementParents[eN_parent] = eN_parent;
6188 for(
int nN_element_0=0,nN_midpoint=0;nN_element_0<multilevelMesh.
meshArray[i-1].
nNodes_element;nN_element_0++)
6189 for(
int nN_element_1=nN_element_0+1;nN_element_1<multilevelMesh.
meshArray[i-1].
nNodes_element;nN_element_1++,nN_midpoint++)
6193 midpoints[nN_midpoint]);
6194 nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
6195 if(nodeItr == newNodeSet.end())
6197 midpoints[nN_midpoint].
nN = nN_new;
6198 newNodeSet.insert(midpoints[nN_midpoint]);
6202 midpoints[nN_midpoint].
nN = nodeItr->nN;
6205 valarray<int> newElement(4);
6206 newElement[0] = eN_new;
6208 newElement[2] = midpoints[0].
nN;
6209 newElement[3] = midpoints[1].
nN;
6210 newElements.push_back(newElement);
6212 newElement[0] = eN_new;
6214 newElement[2] = midpoints[0].
nN;
6215 newElement[3] = midpoints[2].
nN;
6216 newElements.push_back(newElement);
6218 newElement[0] = eN_new;
6220 newElement[2] = midpoints[1].
nN;
6221 newElement[3] = midpoints[2].
nN;
6222 newElements.push_back(newElement);
6224 newElement[0] = eN_parent;
6225 newElement[1] = midpoints[0].
nN;
6226 newElement[2] = midpoints[1].
nN;
6227 newElement[3] = midpoints[2].
nN;
6228 newElements.push_back(newElement);
6236 ebN_neighbor_element=0;
6241 if (eN_neighbor != -1 && elementsForUniform.find(eN_neighbor) == elementsForUniform.end())
6243 if (elementsForTrisection.find(eN_neighbor) == elementsForTrisection.end())
6245 if(elementsForBisection.find(eN_neighbor) == elementsForBisection.end())
6248 int leftNode,rightNode,
6249 leftNode1,rightNode1,
6250 ebN_neighbor_element1,ebN_global1,
6251 longestEdge,longestEdge_element;
6256 longestEdge=ebN_global;
6257 longestEdge_element=ebN_neighbor_element;
6262 ebN_neighbor_element1];
6268 longestEdge = ebN_global1;
6269 longestEdge_element=ebN_neighbor_element1;
6272 if (longestEdge == ebN_global)
6274 elementsForBisection[eN_neighbor] = longestEdge_element;
6278 newElementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6279 newElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6280 elementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6281 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6286 if (elementsForBisection[eN_neighbor] != ebN_neighbor_element)
6288 newElementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6289 newElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6290 elementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6291 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6292 elementsForBisection.erase(eN_neighbor);
6298 if (elementsForTrisection[eN_neighbor][0] != ebN and elementsForTrisection[eN_neighbor][0] != ebN)
6300 nextElementsForUniform.insert(eN_neighbor);
6301 elementsForUniform.insert(eN_neighbor);
6302 elementsForTrisection.erase(eN_neighbor);
6303 newElementsForTrisection.erase(eN_neighbor);
6309 newElementsForUniform.clear();
6310 newElementsForUniform = nextElementsForUniform;
6311 nextElementsForUniform.clear();
6313 for(map<
int,
vector<int> >::iterator eN_trisection_itr = newElementsForTrisection.begin(); eN_trisection_itr != newElementsForTrisection.end(); eN_trisection_itr++)
6315 int eN_parent = eN_trisection_itr->first,
6316 longestEdge_element = eN_trisection_itr->second[0],
6317 otherEdge_element = eN_trisection_itr->second[1];
6325 nodeItr = newNodeSet.find(midpoint_new);
6326 if(nodeItr == newNodeSet.end())
6328 midpoint_new.
nN = nN_new;
6329 newNodeSet.insert(midpoint_new);
6335 longestEdge_element],
6336 ebN_neighbor_element=0;
6341 if (eN_neighbor != -1 && elementsForUniform.find(eN_neighbor) == elementsForUniform.end())
6343 if (elementsForTrisection.find(eN_neighbor) == elementsForTrisection.end())
6345 if(elementsForBisection.find(eN_neighbor) == elementsForBisection.end())
6348 int leftNode,rightNode,
6349 leftNode1,rightNode1,
6350 ebN_neighbor_element1,ebN_global1,
6351 longestEdge,longestEdge_element;
6356 longestEdge=ebN_global;
6357 longestEdge_element=ebN_neighbor_element;
6362 ebN_neighbor_element1];
6369 longestEdge = ebN_global1;
6370 longestEdge_element=ebN_neighbor_element1;
6373 if (longestEdge == ebN_global)
6375 elementsForBisection[eN_neighbor] = longestEdge_element;
6379 nextElementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6380 nextElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6381 elementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6382 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6387 assert(elementsForBisection[eN_neighbor] != ebN_neighbor_element);
6388 nextElementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6389 nextElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6390 elementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6391 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6392 elementsForBisection.erase(eN_neighbor);
6397 if (elementsForTrisection[eN_neighbor][0] != longestEdge and elementsForTrisection[eN_neighbor][0] != longestEdge)
6399 newElementsForUniform.insert(eN_neighbor);
6400 elementsForUniform.insert(eN_neighbor);
6401 elementsForTrisection.erase(eN_neighbor);
6402 nextElementsForTrisection.erase(eN_neighbor);
6407 Node midpoint_other;
6411 nodeItr = newNodeSet.find(midpoint_other);
6412 assert(nodeItr != newNodeSet.end());
6414 newElementsForTrisection.clear();
6415 newElementsForTrisection = nextElementsForTrisection;
6416 nextElementsForTrisection.clear();
6419 for(map<
int,
vector<int> >::iterator eN_trisection_itr = elementsForTrisection.begin(); eN_trisection_itr != elementsForTrisection.end(); eN_trisection_itr++)
6421 int eN_parent = eN_trisection_itr->first,
6422 longestEdge_element = eN_trisection_itr->second[0],
6423 otherEdge_element = eN_trisection_itr->second[1];
6431 if (otherEdge_rightNode == longestEdge_leftNode)
6433 n0 = otherEdge_rightNode;
6434 n1 = otherEdge_leftNode;
6435 n2 = longestEdge_rightNode;
6437 else if (otherEdge_rightNode == longestEdge_rightNode)
6439 n0 = otherEdge_rightNode;
6440 n1 = otherEdge_leftNode;
6441 n2 = longestEdge_leftNode;
6443 else if (otherEdge_leftNode == longestEdge_leftNode)
6445 n0 = otherEdge_leftNode;
6446 n1 = otherEdge_rightNode;
6447 n2 = longestEdge_rightNode;
6451 n0 = otherEdge_leftNode;
6452 n1 = otherEdge_rightNode;
6453 n2 = longestEdge_leftNode;
6455 oldElements.erase(eN_parent);
6456 elementChildren[eN_parent].push_back(eN_new+0);
6457 elementChildren[eN_parent].push_back(eN_new+1);
6458 elementChildren[eN_parent].push_back(eN_parent);
6459 elementParents[eN_new + 0] = eN_parent;
6460 elementParents[eN_new + 1] = eN_parent;
6461 elementParents[eN_parent] = eN_parent;
6467 nodeItr = newNodeSet.find(midpoint_new);
6468 assert(nodeItr != newNodeSet.end());
6469 midpoint_new.
nN = nodeItr->nN;
6470 Node midpoint_other;
6474 nodeItr = newNodeSet.find(midpoint_other);
6475 assert(nodeItr != newNodeSet.end());
6476 midpoint_other.
nN = nodeItr->nN;
6478 valarray<int> newElement(4);
6479 newElement[0] = eN_new;
6481 newElement[2] = midpoint_new.
nN;
6482 newElement[3] = midpoint_other.
nN;
6483 newElements.push_back(newElement);
6485 newElement[0] = eN_new;
6487 newElement[2] = midpoint_other.
nN;
6488 newElement[3] = midpoint_new.
nN;
6489 newElements.push_back(newElement);
6491 newElement[0] = eN_parent;
6493 newElement[2] = midpoint_new.
nN;
6495 newElements.push_back(newElement);
6500 for(map<int,int>::iterator eN_bisect_itr = elementsForBisection.begin(); eN_bisect_itr != elementsForBisection.end(); eN_bisect_itr++)
6502 int eN_parent = eN_bisect_itr->first;
6534 oldElements.erase(eN_parent);
6536 elementChildren[eN_parent].push_back(eN_new);
6537 elementChildren[eN_parent].push_back(eN_parent);
6538 elementParents[eN_new] = eN_parent;
6539 elementParents[eN_parent] = eN_parent;
6542 int nN_element_1,nN_element_2;
6543 bool foundMidpoint=
false;
6546 foundMidpoint=
false;
6547 nN_element_1 = (nN_element_0+1)%3;
6548 nN_element_2 = (nN_element_0+2)%3;
6552 nodeItr = newNodeSet.find(m);
6553 if(nodeItr != newNodeSet.end())
6556 valarray<int> newElement(4);
6557 newElement[0] = eN_new;
6560 newElement[3] = nodeItr->nN;
6561 newElements.push_back(newElement);
6563 newElement[0] = eN_parent;
6565 newElement[2] = nodeItr->nN;
6567 newElements.push_back(newElement);
6585 for(set<int>::iterator eN_itr=oldElements.begin();eN_itr != oldElements.end();eN_itr++)
6593 for(
vector<valarray<int> >::iterator element_itr=newElements.begin();element_itr != newElements.end();element_itr++)
6595 int eN = (*element_itr)[0];
6610 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
6627 for(
unsigned int childE=0;childE<elementChildren[eN_parent].size();childE++)
6647 int * elementTagArray)
6649 using namespace std;
6655 int nLevelsPrev = multilevelMesh.
nLevels;
6657 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
6658 int nElements_tagged = 0;
6661 if (elementTagArray[eN] > 0)
6673 for (
int eN_parent = 0;
6677 if (elementTagArray[eN_parent] == 1)
6703 for (
int eN_parent = 0;
6706 int nBisectedEdges = 0;
6707 if (refined[eN_parent])
6709 for (
int ebN_local = 0;
6716 if (edgeMidNodesArray[ebN] >= 0)
6720 if (nBisectedEdges == 3)
6722 else if (nBisectedEdges == 2)
6724 else if (nBisectedEdges == 1)
6769 if (edgeMidNodesArray[ebN_parent] >= 0)
6775 midpoints[0].
nN = edgeMidNodesArray[ebN_parent];
6787 bool subdivideFailed =
false;
6788 for (
int eN_parent = 0;
6802 if (subdivideFailed)
6825 int& nElements_global,
6827 std::vector<double>& nodeArray,
6828 std::vector<int>& elementNodesArray,
6829 std::vector<int>& elementNeighborsArray,
6830 std::vector<std::list<int> >& childrenList,
6831 std::vector<int>& elementParentsArray,
6832 std::vector<int>& bases,
6833 std::vector<bool>& refined)
6836 bool failed =
false;
6837 const int nElementBoundaries_element = 3;
6838 const int nSpace = 2;
6845 double x[3] = {0.0,0.0,0.0};
6853 int ebN_base = bases[eN];
6854 ib[0] = ebN_base; ib[1] = (ebN_base+1)%nElementBoundaries_element; ib[2]=(ebN_base+2)%nElementBoundaries_element;
6855 for (
int i=0; i < nElementBoundaries_element; i++)
6856 IB[i] = elementNodesArray[nElementBoundaries_element*eN + ib[i]];
6859 int eN_neig = elementNeighborsArray[nElementBoundaries_element*eN + ib[0]];
6860 assert(eN_neig < nElements_global);
6868 int newNodeNumber = nNodes_global;
6869 x[0] = 0.0; x[1] = 0.0; x[2] = 0.0;
6870 for (
int nN = 0; nN < nSpace; nN++)
6871 for (
int I = 0; I < 3; I++)
6872 x[I] += 0.5*nodeArray[3*IB[nN+1]+I];
6879 for (
int I = 0; I < 3; I++)
6880 nodeArray.push_back(x[I]);
6892 int newElementNumber = nElements_global;
6893 E1[0] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
6894 E1[1] = elementNodesArray[nElementBoundaries_element*eN + ib[1]];
6895 E1[2] = newNodeNumber;
6896 E2[0] = newNodeNumber;
6897 E2[1] = elementNodesArray[nElementBoundaries_element*eN + ib[2]];
6898 E2[2] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
6901 N1[1] = newElementNumber;
6902 N1[2] = elementNeighborsArray[nElementBoundaries_element*eN + ib[2]];
6905 N2[0] = elementNeighborsArray[nElementBoundaries_element*eN + ib[1]];
6910 assert(N2[0] < nElements_global);
6911 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
6913 if (elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] == eN)
6915 elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] =
6935 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
6936 elementNodesArray[nElementBoundaries_element*eN + ebN]=E1[ebN];
6938 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
6939 elementNodesArray.push_back(E2[ebN]);
6941 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
6942 elementNeighborsArray[nElementBoundaries_element*eN + ebN]=N1[ebN];
6943 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
6944 elementNeighborsArray.push_back(N2[ebN]);
6951 if (
unsigned(eN) < refined.size() &&
6954 elementParentsArray[eN] = eN;
6955 elementParentsArray.push_back(eN);
6956 childrenList[eN].push_back(eN);
6957 childrenList[eN].push_back(newElementNumber);
6962 assert(
unsigned(eN) < elementParentsArray.size());
6964 elementParentsArray.push_back(elementParentsArray[eN]);
6967 while (
unsigned(eN_tmp) >= refined.size())
6968 eN_tmp = elementParentsArray[eN_tmp];
6969 assert(
unsigned(elementParentsArray[eN_tmp]) < childrenList.size());
6971 childrenList[eN_tmp].push_back(newElementNumber);
6976 nElements_global += 1;
6979 if (
unsigned(eN) < refined.size())
6982 else if (elementNeighborsArray[nElementBoundaries_element*eN_neig + bases[eN_neig]] == eN)
6990 int newNodeNumber = nNodes_global;
6991 x[0] = 0.0; x[1] = 0.0; x[2] = 0.0;
6992 for (
int nN = 0; nN < nSpace; nN++)
6993 for (
int I = 0; I < 3; I++)
6995 x[I] += 0.5*nodeArray[3*IB[nN+1]+I];
6998 for (
int I = 0; I < 3; I++)
6999 nodeArray.push_back(x[I]);
7018 int ebN_base_neig = bases[eN_neig];
7019 ibn[0]=ebN_base_neig; ibn[1]=(ebN_base_neig+1)%nElementBoundaries_element;
7020 ibn[2]=(ebN_base_neig+2)%nElementBoundaries_element;
7022 int newElementNumber = nElements_global;
7023 int newElementNumberNeig = nElements_global+1;
7024 E1[0] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
7025 E1[1] = elementNodesArray[nElementBoundaries_element*eN + ib[1]];
7026 E1[2] = newNodeNumber;
7027 E2[0] = newNodeNumber;
7028 E2[1] = elementNodesArray[nElementBoundaries_element*eN + ib[2]];
7029 E2[2] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
7031 N1[0] = newElementNumberNeig;
7032 N1[1] = newElementNumber;
7033 N1[2] = elementNeighborsArray[nElementBoundaries_element*eN + ib[2]];
7037 N2[0] = elementNeighborsArray[nElementBoundaries_element*eN + ib[1]];
7042 assert(N2[0] < nElements_global);
7043 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7045 if (elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] == eN)
7047 elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] =
7056 E1n[0] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[0]];
7057 E1n[1] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[1]];
7058 E1n[2] = newNodeNumber;
7059 E2n[0] = newNodeNumber;
7060 E2n[1] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[2]];
7061 E2n[2] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[0]];
7063 N1n[0] = newElementNumber;
7064 N1n[1] = newElementNumberNeig;
7065 N1n[2] = elementNeighborsArray[nElementBoundaries_element*eN_neig + ibn[2]];
7069 N2n[0] = elementNeighborsArray[nElementBoundaries_element*eN_neig + ibn[1]];
7074 assert(N2n[0] < nElements_global);
7075 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7077 if (elementNeighborsArray[nElementBoundaries_element*N2n[0] + ebN] == eN_neig)
7079 elementNeighborsArray[nElementBoundaries_element*N2n[0] + ebN] =
7080 newElementNumberNeig;
7088 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7090 elementNodesArray[nElementBoundaries_element*eN + ebN] =E1[ebN];
7091 elementNodesArray[nElementBoundaries_element*eN_neig + ebN]=E1n[ebN];
7094 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7095 elementNodesArray.push_back(E2[ebN]);
7097 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7098 elementNodesArray.push_back(E2n[ebN]);
7101 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7103 elementNeighborsArray[nElementBoundaries_element*eN + ebN] =N1[ebN];
7104 elementNeighborsArray[nElementBoundaries_element*eN_neig + ebN]=N1n[ebN];
7106 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7107 elementNeighborsArray.push_back(N2[ebN]);
7108 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7109 elementNeighborsArray.push_back(N2n[ebN]);
7118 if (
unsigned(eN) < refined.size() &&
7121 elementParentsArray[eN] = eN;
7122 elementParentsArray.push_back(eN);
7123 childrenList[eN].push_back(eN);
7124 childrenList[eN].push_back(newElementNumber);
7129 assert(
unsigned(eN) < elementParentsArray.size());
7131 elementParentsArray.push_back(elementParentsArray[eN]);
7134 while (
unsigned(eN_tmp) >= refined.size())
7135 eN_tmp = elementParentsArray[eN_tmp];
7136 assert(
unsigned(elementParentsArray[eN_tmp]) < childrenList.size());
7138 childrenList[eN_tmp].push_back(newElementNumber);
7141 if (
unsigned(eN_neig) < refined.size() &&
7144 elementParentsArray[eN_neig] = eN_neig;
7145 elementParentsArray.push_back(eN_neig);
7146 childrenList[eN_neig].push_back(eN_neig);
7147 childrenList[eN_neig].push_back(newElementNumberNeig);
7152 assert(
unsigned(eN_neig) < elementParentsArray.size());
7154 elementParentsArray.push_back(elementParentsArray[eN_neig]);
7156 int eN_tmp = eN_neig;
7157 while (
unsigned(eN_tmp) >= refined.size())
7158 eN_tmp = elementParentsArray[eN_tmp];
7159 assert(
unsigned(elementParentsArray[eN_tmp]) < childrenList.size());
7161 childrenList[eN_tmp].push_back(newElementNumberNeig);
7165 nElements_global += 2;
7168 if (
unsigned(eN) < refined.size())
7171 if (
unsigned(eN_neig) < refined.size())
7172 refined[eN_neig] =
true;
7183 elementNeighborsArray,
7185 elementParentsArray,
7193 elementNeighborsArray,
7195 elementParentsArray,
7204 std::vector<bool>& refined,
7205 std::vector<int>& edgeMidNodesArray,
7206 const int* elementNodesArray,
7207 const int* elementBoundariesArray,
7208 const int* elementNeighborsArray,
7209 const double * nodeArray)
7211 bool failed =
false;
7212 const int nElementBoundaries_element = 3;
7213 bool refinedAlready = refined[eN];
7222 eN_longest = elementBoundariesArray[eN*nElementBoundaries_element+eN_longest_local];
7223 for (
int ebN_local = 0; ebN_local < nElementBoundaries_element; ebN_local++)
7225 const int ebN = elementBoundariesArray[eN*nElementBoundaries_element+ebN_local];
7227 if (edgeMidNodesArray[ebN] >= 0 && !refinedAlready)
7234 assert(edgeMidNodesArray[ebN] < 0 || refinedAlready);
7235 if (edgeMidNodesArray[ebN] < 0)
7237 edgeMidNodesArray[ebN] = nNodes_global;
7240 int eN_neig = elementNeighborsArray[eN*nElementBoundaries_element+ebN_local];
7250 elementBoundariesArray,
7251 elementNeighborsArray,
7261 int ebN_neig,
int eN_neig,
7263 std::vector<bool>& refined,
7264 std::vector<int>& edgeMidNodesArray,
7265 const int* elementNodesArray,
7266 const int* elementBoundariesArray,
7267 const int* elementNeighborsArray,
7268 const double * nodeArray)
7274 bool failed =
false;
7276 const int nElementBoundaries_element = 3;
7293 if (edgeMidNodesArray[ebN_neig] < 0)
7295 edgeMidNodesArray[ebN_neig] = nNodes_global;
7300 assert(eN_neig >=0);
7304 int ebN_neig_longest = elementBoundariesArray[eN_neig*nElementBoundaries_element+eN_neig_longest_local];
7306 refined[eN_neig] =
true;
7307 if (edgeMidNodesArray[ebN_neig_longest] < 0)
7309 edgeMidNodesArray[ebN_neig_longest] = nNodes_global;
7312 if (ebN_longest == ebN_neig_longest)
7315 assert(ebN_longest == ebN_neig);
7319 int eN_neig_neig = elementNeighborsArray[eN_neig*nElementBoundaries_element+eN_neig_longest_local];
7321 ebN_neig_longest,eN_neig_neig,
7326 elementBoundariesArray,
7327 elementNeighborsArray,
7332 const int* elementNodesArray,
7333 const double * nodeArray)
7335 const int nElementBoundaries_element = 3;
7336 int longest = 0;
double h_longest=0.0;
7338 int nN0 = elementNodesArray[eN*nElementBoundaries_element+((ebN+1)%nElementBoundaries_element)];
7339 int nN1 = elementNodesArray[eN*nElementBoundaries_element+((ebN+2)%nElementBoundaries_element)];
7340 double len = fabs((nodeArray[3*nN1+0]-nodeArray[3*nN0+0])*(nodeArray[3*nN1+0]-nodeArray[3*nN0+0])+
7341 (nodeArray[3*nN1+1]-nodeArray[3*nN0+1])*(nodeArray[3*nN1+1]-nodeArray[3*nN0+1])+
7342 (nodeArray[3*nN1+2]-nodeArray[3*nN0+2])*(nodeArray[3*nN1+2]-nodeArray[3*nN0+2]));
7344 longest = ebN; h_longest = len;
7345 for (ebN = 1; ebN < nElementBoundaries_element; ebN++)
7347 nN0 = elementNodesArray[eN*nElementBoundaries_element+((ebN+1)%nElementBoundaries_element)];
7348 nN1 = elementNodesArray[eN*nElementBoundaries_element+((ebN+2)%nElementBoundaries_element)];
7349 len = fabs((nodeArray[3*nN1+0]-nodeArray[3*nN0+0])*(nodeArray[3*nN1+0]-nodeArray[3*nN0+0])+
7350 (nodeArray[3*nN1+1]-nodeArray[3*nN0+1])*(nodeArray[3*nN1+1]-nodeArray[3*nN0+1])+
7351 (nodeArray[3*nN1+2]-nodeArray[3*nN0+2])*(nodeArray[3*nN1+2]-nodeArray[3*nN0+2]));
7353 if (len > h_longest)
7354 { len = h_longest; longest = ebN;}
7361 int* elementParentsArray,
7362 int* elementChildrenOffsets,
7363 int* elementChildrenArray,
7364 int* elementNodesArray_child,
7365 const std::vector<int>& edgeMidNodesArray,
7366 const std::vector<bool>& refined,
7367 const int* elementNodesArray_parent,
7368 const int* elementBoundariesArray_parent,
7369 const double* nodeArray_parent)
7371 bool failed =
false,u4T=
false;
7373 const int simplexDim = 3;
7374 const int childOffset = elementChildrenOffsets[eN_parent];
7375 if (!refined[eN_parent])
7378 for (
int nN = 0; nN < simplexDim; nN++)
7380 elementNodesArray_child[eN_parent*simplexDim+nN] =
7381 elementNodesArray_parent[eN_parent*simplexDim+nN];
7384 elementParentsArray[eN_parent] = eN_parent;
7385 elementChildrenOffsets[eN_parent+1]= childOffset+1;
7386 elementChildrenArray[childOffset] = eN_parent;
7391 int nBisectedEdges = 0;
7392 int midnodes[3],vnodes[3];
7393 for (
int ebN_local = 0; ebN_local < simplexDim; ebN_local++)
7395 const int ebN = elementBoundariesArray_parent[eN_parent*simplexDim+ebN_local];
7396 vnodes[ebN_local] = elementNodesArray_parent[eN_parent*simplexDim+ebN_local];
7397 midnodes[ebN_local] = edgeMidNodesArray[ebN];
7398 if (edgeMidNodesArray[ebN] >= 0)
7401 assert(nBisectedEdges > 0);
7403 if (nBisectedEdges == 3)
7410 elementChildrenOffsets[eN_parent+1]= childOffset+4;
7414 elementNodesArray_child[eN_parent*simplexDim+0]=
7416 elementNodesArray_child[eN_parent*simplexDim+1]=
7418 elementNodesArray_child[eN_parent*simplexDim+2]=
7421 elementParentsArray[eN_parent] = eN_parent;
7422 elementChildrenArray[childOffset+0] = eN_parent;
7423 for (
int c=0;
c<simplexDim;
c++)
7425 elementNodesArray_child[eN_new*simplexDim+0]=
7427 elementNodesArray_child[eN_new*simplexDim+1]=
7428 midnodes[(
c+1)%simplexDim];
7429 elementNodesArray_child[eN_new*simplexDim+2]=
7430 midnodes[(
c+2)%simplexDim];
7432 elementParentsArray[eN_new] = eN_parent;
7433 elementChildrenArray[childOffset+
c+1] = eN_new;
7453 elementNodesArray_parent,
7455 assert(midnodes[base] >= 0);
7457 elementChildrenOffsets[eN_parent+1]= childOffset+4;
7460 elementNodesArray_child[eN_parent*simplexDim+0]=
7462 elementNodesArray_child[eN_parent*simplexDim+1]=
7463 midnodes[(base+2)%simplexDim];
7464 elementNodesArray_child[eN_parent*simplexDim+2]=
7467 elementParentsArray[eN_parent] = eN_parent;
7468 elementChildrenArray[childOffset+0] = eN_parent;
7470 elementNodesArray_child[eN_new*simplexDim+0]=
7471 midnodes[(base+2)%simplexDim];
7472 elementNodesArray_child[eN_new*simplexDim+1]=
7473 vnodes[(base+1)%simplexDim];
7474 elementNodesArray_child[eN_new*simplexDim+2]=
7477 elementParentsArray[eN_new] = eN_parent;
7478 elementChildrenArray[childOffset+1] = eN_new;
7481 elementNodesArray_child[eN_new*simplexDim+0]=
7483 elementNodesArray_child[eN_new*simplexDim+1]=
7484 vnodes[(base+2)%simplexDim];
7485 elementNodesArray_child[eN_new*simplexDim+2]=
7486 midnodes[(base+1)%simplexDim];
7488 elementParentsArray[eN_new] = eN_parent;
7489 elementChildrenArray[childOffset+2] = eN_new;
7492 elementNodesArray_child[eN_new*simplexDim+0]=
7493 midnodes[(base+1)%simplexDim];
7494 elementNodesArray_child[eN_new*simplexDim+1]=
7496 elementNodesArray_child[eN_new*simplexDim+2]=
7499 elementParentsArray[eN_new] = eN_parent;
7500 elementChildrenArray[childOffset+3] = eN_new;
7504 else if (nBisectedEdges == 2)
7525 elementNodesArray_parent,
7527 assert(midnodes[base] >= 0);
7529 elementChildrenOffsets[eN_parent+1]= childOffset+3;
7531 if (midnodes[(base+1)%simplexDim] < 0)
7533 assert(midnodes[(base+2)%simplexDim] >= 0);
7536 elementNodesArray_child[eN_parent*simplexDim+0]=
7537 vnodes[(base+2)%simplexDim];
7538 elementNodesArray_child[eN_parent*simplexDim+1]=
7540 elementNodesArray_child[eN_parent*simplexDim+2]=
7543 elementParentsArray[eN_parent] = eN_parent;
7544 elementChildrenArray[childOffset+0] = eN_parent;
7546 elementNodesArray_child[eN_new*simplexDim+0]=
7548 elementNodesArray_child[eN_new*simplexDim+1]=
7549 midnodes[(base+2)%simplexDim];
7550 elementNodesArray_child[eN_new*simplexDim+2]=
7553 elementParentsArray[eN_new] = eN_parent;
7554 elementChildrenArray[childOffset+1] = eN_new;
7557 elementNodesArray_child[eN_new*simplexDim+0]=
7558 midnodes[(base+2)%simplexDim];
7559 elementNodesArray_child[eN_new*simplexDim+1]=
7560 vnodes[(base+1)%simplexDim];
7561 elementNodesArray_child[eN_new*simplexDim+2]=
7564 elementParentsArray[eN_new] = eN_parent;
7565 elementChildrenArray[childOffset+2] = eN_new;
7570 assert(midnodes[(base+2)%simplexDim] < 0);
7573 elementNodesArray_child[eN_parent*simplexDim+0]=
7575 elementNodesArray_child[eN_parent*simplexDim+1]=
7576 vnodes[(base+1)%simplexDim];
7577 elementNodesArray_child[eN_parent*simplexDim+2]=
7580 elementParentsArray[eN_parent] = eN_parent;
7581 elementChildrenArray[childOffset+0] = eN_parent;
7583 elementNodesArray_child[eN_new*simplexDim+0]=
7585 elementNodesArray_child[eN_new*simplexDim+1]=
7587 elementNodesArray_child[eN_new*simplexDim+2]=
7588 midnodes[(base+1)%simplexDim];
7590 elementParentsArray[eN_new] = eN_parent;
7591 elementChildrenArray[childOffset+1] = eN_new;
7594 elementNodesArray_child[eN_new*simplexDim+0]=
7596 elementNodesArray_child[eN_new*simplexDim+1]=
7597 vnodes[(base+2)%simplexDim];
7598 elementNodesArray_child[eN_new*simplexDim+2]=
7599 midnodes[(base+1)%simplexDim];
7601 elementParentsArray[eN_new] = eN_parent;
7602 elementChildrenArray[childOffset+2] = eN_new;
7606 else if (nBisectedEdges == 1)
7618 elementNodesArray_parent,
7621 elementChildrenOffsets[eN_parent+1]= childOffset+2;
7623 assert(midnodes[base] >= 0);
7625 elementNodesArray_child[eN_parent*simplexDim+0]=
7627 elementNodesArray_child[eN_parent*simplexDim+1]=
7628 vnodes[(base+1)%simplexDim];
7629 elementNodesArray_child[eN_parent*simplexDim+2]=
7632 elementParentsArray[eN_parent] = eN_parent;
7633 elementChildrenArray[childOffset+0] = eN_parent;
7635 elementNodesArray_child[eN_new*simplexDim+0]=
7637 elementNodesArray_child[eN_new*simplexDim+1]=
7638 vnodes[(base+2)%simplexDim];
7639 elementNodesArray_child[eN_new*simplexDim+2]=
7642 elementParentsArray[eN_new] = eN_parent;
7643 elementChildrenArray[childOffset+1] = eN_new;
7648 int child_start = elementChildrenOffsets[eN_parent];
7649 int child_end = elementChildrenOffsets[eN_parent+1];
7665#include <pybind11/pybind11.h>
7666#include <pybind11/numpy.h>
7667namespace py = pybind11;
7675 py::module_ meshTools = py::module_::import(
"proteus.MeshTools");
7676 py::object plexMesh = py::reinterpret_borrow<py::object>(dmplexMesh);
7678 if (py::isinstance(plexMesh, meshTools.attr(
"Mesh"))) {
7679 logEvent(
"DMPlex mesh is available in MeshTools.", 1);
7680 logEvent(
"Reading DMPlex mesh from MeshTools", 1);
7682 mesh.
nNodes_global = plexMesh.attr(
"nNodes_global").cast<
int>();
7685 mesh.
nEdges_global = plexMesh.attr(
"nEdges_global").cast<
int>();
7686 mesh.
nNodes_element = plexMesh.attr(
"nNodes_element").cast<
int>();
7694 py::array_t<int> nodeElementsArray = plexMesh.attr(
"nodeElementsArray").cast<py::array_t<int>>();
7695 auto nodeElementsArrayBuf = nodeElementsArray.request();
7698 py::array_t<int> elementNodesArray = plexMesh.attr(
"elementNodesArray").cast<py::array_t<int>>();
7699 auto elementNodesArrayBuf = elementNodesArray.request();
7702 py::array_t<int> nodeElementOffsets = plexMesh.attr(
"nodeElementOffsets").cast<py::array_t<int>>();
7703 auto nodeElementOffsetsBuf = nodeElementOffsets.request();
7710 py::array_t<int> elementBoundariesArray = plexMesh.attr(
"elementBoundariesArray").cast<py::array_t<int>>();
7711 auto elementBoundariesArrayBuf = elementBoundariesArray.request();
7725 py::array_t<int> elementBoundaryElementsArray = plexMesh.attr(
"elementBoundaryElementsArray").cast<py::array_t<int>>();
7726 auto elementBoundaryElementsArrayBuf = elementBoundaryElementsArray.request();
7733 py::array_t<int> interiorElementBoundariesArray = plexMesh.attr(
"interiorElementBoundariesArray").cast<py::array_t<int>>();
7734 auto interiorElementBoundariesArrayBuf = interiorElementBoundariesArray.request();
7737 py::array_t<int> exteriorElementBoundariesArray = plexMesh.attr(
"exteriorElementBoundariesArray").cast<py::array_t<int>>();
7738 auto exteriorElementBoundariesArrayBuf = exteriorElementBoundariesArray.request();
7741 py::array_t<int> edgeNodesArray = plexMesh.attr(
"edgeNodesArray").cast<py::array_t<int>>();
7742 auto edgeNodesArrayBuf = edgeNodesArray.request();
7745 py::array_t<int> nodeStarArray = plexMesh.attr(
"nodeStarArray").cast<py::array_t<int>>();
7746 auto nodeStarArrayBuf = nodeStarArray.request();
7747 mesh.
nodeStarArray =
static_cast<int*
>(nodeStarArrayBuf.ptr);
7749 py::array_t<int> nodeStarOffsets = plexMesh.attr(
"nodeStarOffsets").cast<py::array_t<int>>();
7750 auto nodeStarOffsetsBuf = nodeStarOffsets.request();
7753 py::array_t<int> elementMaterialTypes = plexMesh.attr(
"elementMaterialTypes").cast<py::array_t<int>>();
7754 auto elementMaterialTypesBuf = elementMaterialTypes.request();
7757 py::array_t<int> elementBoundaryMaterialTypes = plexMesh.attr(
"elementBoundaryMaterialTypes").cast<py::array_t<int>>();
7758 auto elementBoundaryMaterialTypesBuf = elementBoundaryMaterialTypes.request();
7761 py::array_t<int> nodeMaterialTypes = plexMesh.attr(
"nodeMaterialTypes").cast<py::array_t<int>>();
7762 auto nodeMaterialTypesBuf = nodeMaterialTypes.request();
7765 py::array_t<double> nodeArray = plexMesh.attr(
"nodeArray").cast<py::array_t<double>>();
7766 auto nodeArrayBuf = nodeArray.request();
7767 mesh.
nodeArray =
static_cast<double*
>(nodeArrayBuf.ptr);
7803 std::cerr <<
"DMPlex mesh is not available in MeshTools." << std::endl;
7809#ifdef MWF_HACK_2D_COARSE
7810int findGlueNeighbor2d(
int eN,
7811 std::vector<int>& elementNeighborsArray,
7812 std::vector<int>& elementParentsArray)
7817 const int nElementBoundaries_element = 3;
7819 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7821 const int eN_neig = elementNeighborsArray[eN*nElementBoundaries_element+ebN];
7822 if (eN_neig >= 0 && elementParentsArray[eN_neig] == elementParentsArray[eN])
7830int findT2Neighbor(
int eN,
int eN_base,
7831 std::vector<int>& elementNeighborsArray,
7832 std::vector<int>& elementParentsArray)
7841 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7843 const int eN_neig = elementNeighborsArray[eN*nElementBoundaries_element+ebN];
7844 if (ebN != eN_base && eN_neig >= 0 &&
7845 elementParentsArray[eN_neig] != elementParentsArray[eN])
7853bool newestNodeGlue(
int eN,
7854 std::vector<bool>& mayCoarsen,
7855 int& nElements_global,
7857 std::vector<double>& nodeArray,
7858 std::vector<int>& elementNodesArray,
7859 std::vector<int>& elementNeighborsArray,
7860 std::vector<int>& elementParentsArray,
7861 std::vector<int>& bases,
7862 std::vector<bool>& coarsened)
7864 using namespace std;
7865 bool failed =
false;
7866 if (!mayCoarsen[eN])
7867 {failed =
true;
return failed;}
7869 const int nElementBoundaries_element = 3;
7870 bool connectable =
false;
7871 int eN_base = bases[eN];
7872 int nN_global_base = elementNodesArray[eN*nElementBoundaries_element+eN_base];
7873 int eN_glue = -1, eN_t2 = -1;
7874 eN_glue = findGlueNeighbor(eN,
7875 elementNeighborsArray,
7876 elementParentsArray);
7877 eN_t2 = findT2Neighbor(eN,eN_base,
7878 elementNeighborsArray,
7879 elementParentsArray);
7881 {failed =
true;
return failed;}
7882 if (eN_t2 >= 0 && elementNodesArray[eN_t2*nElementBoundaries_element+bases[eN_t2]] != nN_global_base)
7883 failed = newestNodeGlue(eN_t2,
7889 elementNeighborsArray,
7890 elementParentsArray,
7896 connectable = (nN_global_base == elementNodesArray[eN_glue*nElementBoundaries_element+bases[eN_glue]]);
7898 failed = newestNodeGlue(eN_glue,
7904 elementNeighborsArray,
7905 elementParentsArray,
7911 eN_glue = findGlueNeighbor(eN,
7912 elementNeighborsArray,
7913 elementParentsArray);
7915 { failed =
true;
return failed;}
7916 eN_t2 = findT2Neighbor(eN,eN_base,
7917 elementNeighborsArray,
7918 elementParentsArray);
7919 bool t2connectable =
false;
7922 int eN_t2_glue = -1;
7924 eN_t2_glue = findGlueNeighbor(eN_t2,
7925 elementNeighborsArray,
7926 elementParentsArray);
7928 if (eN_t2_glue < 0 && eN_t2 >= 0)
7929 {failed =
true;
return failed;}
7930 t2connectable = (eN_t2 >= 0 &&
7931 (elementNodesArray[eN_t2*nElementBoundaries_element+bases[eN_t2]] ==
7932 elementNodesArray[eN_t2_glue*nElementBoundaries_element+bases[eN_t2_glue]]));
7933 if (eN_t2 >= 0 && !t2connectable)
7934 failed = newestNodeGlue(eN_t2_glue,
7940 elementNeighborsArray,
7941 elementParentsArray,
7946 t2connectable =
false;
7948 eN_t2_glue = findGlueNeighbor(eN_t2,
7949 elementNeighborsArray,
7950 elementParentsArray);
7951 if (eN_t2 >= 0 && eN_t2_glue >= 0)
7952 t2connectable = (elementNodesArray[eN_t2*nElementBoundaries_element+bases[eN_t2]] ==
7953 elementNodesArray[eN_t2_glue*nElementBoundaries_element+bases[eN_t2_glue]]);
7955 failed = glueElements(eN_t2,eN_t2_glue,
7961 elementNeighborsArray,
7962 elementParentsArray,
7967 failed = glueElements(eN,eN_glue,
7973 elementNeighborsArray,
7974 elementParentsArray,
7981bool glueElements(
int eN,
int eN_glue,
7982 std::vector<bool>& mayCoarsen,
7983 int& nElements_global,
7985 std::vector<double>& nodeArray,
7986 std::vector<int>& elementNodesArray,
7987 std::vector<int>& elementNeighborsArray,
7988 std::vector<int>& elementParentsArray,
7989 std::vector<int>& bases,
7990 std::vector<bool>& coarsened)
7994 const int nElementBoundaries_element = 3;
7995 assert(eN >= 0 && eN_glue >= 0);
7996 if (!mayCoarsen[eN] || !mayCoarsen[eN_glue])
7998 bool connectable = (elementNodesArray[eN*nElementBoundaries_element+bases[eN]] ==
7999 elementNodesArray[eN_glue*nElementBoundaries_element+bases[eN_glue]]);
8014 base_parent = (bases[eN]+1) % nElementBoundaries_element;
8015 int nN_global_gone = elementNodesArray[eN*nElementBoundaries_element+bases[eN]];
8017 for (
int nN = 0; nN < nElementBoundaries_element; nN++)
8019 if (nN == bases[eN])
8020 elementNodesArray[eN*nElementBoundaries_element+ebN] = (bases[eN_glue]+1) % nElementBoundaries_element;
8022 coarsened[eN] =
true;
float * vector(long nl, long nh)
#define DEFAULT_ELEMENT_MATERIAL
int regularHexahedralToTetrahedralElementBoundaryMaterials(const double &Lx, const double &Ly, const double &Lz, Mesh &mesh)
int newTetrahedron(int eN, int *nodes, int n0, int n1, int n2, int n3)
int writeTriangleMesh(Mesh &mesh, const char *filebase, int triangleIndexBase)
int writeNodes(std::ostream &meshFile, const Mesh &mesh)
int constructElementBoundaryElementsArray_edge(Mesh &mesh)
int allocateGeometricInfo_hexahedron(Mesh &mesh)
int computeGeometricInfo_NURBS(Mesh &mesh)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_hexahedron(Mesh &mesh)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_quadrilateral(Mesh &mesh)
int constructElementBoundaryElementsArray_quadrilateral(Mesh &mesh)
int setNewestNodeBasesToLongestEdge(MultilevelMesh &multilevelMesh)
int allocateGeometricInfo_triangle(Mesh &mesh)
int readTetgenMesh(Mesh &mesh, const char *filebase, int tetgenIndexBase)
int read3DM(Mesh &mesh, const char *filebase, int indexBase)
int writeElements(std::ostream &meshFile, const Mesh &mesh)
int constructElementBoundaryElementsArray_NURBS(Mesh &mesh)
int regularEdgeMeshNodes(const int &nx, const double &Lx, Mesh &mesh)
int readTriangleElementBoundaryMaterialTypes(Mesh &mesh, const char *filebase, int triangleIndexBase)
int locallyRefineTriangleMesh(MultilevelMesh &multilevelMesh, int *elementTagArray)
int regularHexahedralMeshElementBoundaryMaterials(const double &Lx, const double &Ly, const double &Lz, Mesh &mesh)
int regularRectangularToTriangularMeshElements(const int &nx, const int &ny, Mesh &mesh, int triangleFlag)
int readDMPlexMesh(PyObject *dmplexMesh, Mesh &mesh)
int write2dmMesh(Mesh &mesh, const char *filebase, int adhIndexBase)
int newEdge(int eN, int *nodes, int n0, int n1)
int regularMeshNodes2D(const int &nx, const int &ny, const double &Lx, const double &Ly, Mesh &mesh)
int constructElementBoundaryElementsArray_hexahedron(Mesh &mesh)
int newQuadrilateral(int eN, int *nodes, int n0, int n1, int n2, int n3)
void initializeMesh(Mesh &mesh)
int regularQuadrilateralMeshElements(const int &nx, const int &ny, Mesh &mesh)
int read2DM(Mesh &mesh, const char *filebase, int indexBase)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_tetrahedron(Mesh &mesh)
int findLocalLongestEdge2d(int eN, const int *elementNodesArray, const double *nodeArray)
bool add4TnodesForConformity2d(int eN, int ebN_longest, int ebN_neig, int eN_neig, int &nNodes_global, std::vector< bool > &refined, std::vector< int > &edgeMidNodesArray, const int *elementNodesArray, const int *elementBoundariesArray, const int *elementNeighborsArray, const double *nodeArray)
int regularHexahedralToTetrahedralMeshNodes(const int &nx, const int &ny, const int &nz, const double &Lx, const double &Ly, const double &Lz, Mesh &mesh)
int globallyRefineHexahedralMesh(const int &nLevels, Mesh &mesh, MultilevelMesh &multilevelMesh, bool averageNewNodeFlags)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_edge(Mesh &mesh)
int regularMeshNodes(const int &nx, const int &ny, const int &nz, const double &Lx, const double &Ly, const double &Lz, Mesh &mesh)
int readBC(Mesh &mesh, const char *filebase, int indexBase)
int allocateGeometricInfo_quadrilateral(Mesh &mesh)
int constructElementBoundaryElementsArray_tetrahedron(Mesh &mesh)
int regularHexahedralToTetrahedralMeshElements(const int &nx, const int &ny, const int &nz, Mesh &mesh)
int computeGeometricInfo_edge(Mesh &mesh)
int allocateGeometricInfo_edge(Mesh &mesh)
bool add4TnodesForRefinement2d(int eN, int &nNodes_global, std::vector< bool > &refined, std::vector< int > &edgeMidNodesArray, const int *elementNodesArray, const int *elementBoundariesArray, const int *elementNeighborsArray, const double *nodeArray)
int assignElementBoundaryMaterialTypesFromParent(Mesh &parentMesh, Mesh &childMesh, const int *levelElementParentsArray, const int &nSpace_global)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_triangle(Mesh &mesh)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_NURBS(Mesh &mesh)
int newHexahedron(int eN, int *nodes, int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
int newTriangle(int eN, int *nodes, int n0, int n1, int n2)
int constructElementBoundaryElementsArray_triangle(Mesh &mesh)
int regularRectangularToTriangularMeshNodes(const int &nx, const int &ny, const double &Lx, const double &Ly, Mesh &mesh)
int locallyRefineTriangleMesh_redGreen(MultilevelMesh &multilevelMesh, int *elementTagArray)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_quadrilateral(Mesh &mesh)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_triangle(Mesh &mesh)
bool newestNodeBisect(int eN, int &nElements_global, int &nNodes_global, std::vector< double > &nodeArray, std::vector< int > &elementNodesArray, std::vector< int > &elementNeighborsArray, std::vector< std::list< int > > &childrenList, std::vector< int > &elementParentsArray, std::vector< int > &bases, std::vector< bool > &refined)
int computeGeometricInfo_quadrilateral(Mesh &mesh)
int computeGeometricInfo_triangle(Mesh &mesh)
int regularQuadrilateralMeshElementBoundaryMaterials(const double &Lx, const double &Ly, Mesh &mesh)
int readTriangleMesh(Mesh &mesh, const char *filebase, int triangleIndexBase)
int writeTetgenMesh(Mesh &mesh, const char *filebase, int tetgenIndexBase)
int globallyRefineTetrahedralMesh(const int &nLevels, Mesh &mesh, MultilevelMesh &multilevelMesh, bool averageNewNodeFlags)
void midpoint(const double *left, const double *right, Node &midpoint)
int edgeMeshElements(const int &nx, Mesh &mesh)
int globallyRefineEdgeMesh(const int &nLevels, Mesh &mesh, MultilevelMesh &multilevelMesh, bool averageNewNodeFlags)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryNumbers_edge(Mesh &mesh)
bool subdivideTriangle4T(int eN_parent, int &eN_new, int *elementParentsArray, int *elementChildrenOffsets, int *elementChildrenArray, int *elementNodesArray_child, const std::vector< int > &edgeMidNodesArray, const std::vector< bool > &refined, const int *elementNodesArray_parent, const int *elementBoundariesArray_parent, const double *nodeArray_parent)
int allocateNodeAndElementNodeDataStructures(Mesh &mesh, int nElements_global, int nNodes_global, int nNodes_element)
int globallyRefineTriangularMesh(const int &nLevels, Mesh &mesh, MultilevelMesh &multilevelMesh, bool averageNewNodeFlags)
int readElements(std::istream &meshFile, Mesh &mesh)
int readHex(Mesh &mesh, const char *filebase, int indexBase)
double edgeLengthFromNodeNumbers(double *nodeArray, const int &left, const int &right)
int allocateGeometricInfo_NURBS(Mesh &mesh)
int globallyRefineQuadrilateralMesh(const int &nLevels, Mesh &mesh, MultilevelMesh &multilevelMesh, bool averageNewNodeFlags)
int write3dmMesh(Mesh &mesh, const char *filebase, int adhIndexBase)
int allocateGeometricInfo_tetrahedron(Mesh &mesh)
int regularHexahedralMeshElements(const int &nx, const int &ny, const int &nz, const int &px, const int &py, const int &pz, Mesh &mesh)
int regularNURBSMeshElements(const int &nx, const int &ny, const int &nz, const int &px, const int &py, const int &pz, Mesh &mesh)
int regularRectangularToTriangularElementBoundaryMaterials(const double &Lx, const double &Ly, Mesh &mesh)
int locallyRefineTriangleMesh_4T(MultilevelMesh &multilevelMesh, int *elementTagArray)
int constructElementBoundaryElementsArrayWithGivenElementBoundaryAndEdgeNumbers_tetrahedron(Mesh &mesh)
int locallyRefineEdgeMesh(MultilevelMesh &multilevelMesh, int *elementTagArray)
int readTetgenElementBoundaryMaterialTypes(Mesh &mesh, const char *filebase, int tetgenIndexBase)
int reorientTetrahedralMesh(Mesh &mesh)
int computeGeometricInfo_hexahedron(Mesh &mesh)
int computeGeometricInfo_tetrahedron(Mesh &mesh)
void reorientNodes_tet(double *nodeArray, int *nodes)
double edgeLength(int nL, int nR, const double *nodeArray)
const int INTERIOR_ELEMENT_BOUNDARY_MATERIAL
int growMultilevelMesh(int nLevels2add, MultilevelMesh &multilevelMesh)
double tetrahedronVolume(int n0, int n1, int n2, int n3, const double *nodeArray)
const int EXTERIOR_NODE_MATERIAL
const int DEFAULT_NODE_MATERIAL
const int INTERIOR_NODE_MATERIAL
double hexahedronVolume(int n0, int n1, int n2, int n3, int n4, int n5, int n6, int n7, const double *nodeArray)
const int EXTERIOR_ELEMENT_BOUNDARY_MATERIAL
double triangleArea(int n0, int n1, int n2, const double *nodeArray)
bool write3dmMeshNodesAndElements(const char *filebase, const int &indexBase, const int &nElements, const int &nNodes, const double *nodeArray, const int *elementNodesArray, const int *elementMaterialTypes)
bool writeTriangleElementBoundaryNodes(const char *filebase, const int &indexBase, const int &nElementBoundaries, const int *elementBoundaryNodesArray, const int *elementBoundaryMaterialTypes)
bool readTetgenElementBoundaries(const char *filebase, const int &indexBase, bool &hasMarkers, int &nElementBoundaries, std::vector< int > &elementBoundaryNodesArray, std::vector< int > &elementBoundaryMaterialTypesArray, const int &defaultBoundaryMaterialType)
bool readTriangleMeshNodesAndElements(const char *filebase, const int &indexBase, int &nElements, int &nNodes, std::vector< double > &nodeArray, std::vector< int > &elementNodesArray, std::vector< int > &nodeMaterialTypes, std::vector< int > &elementMaterialTypes, const int &defaultElementMaterialType, const int &defaultNodeMaterialType)
bool writeTriangleMeshNodesAndElements(const char *filebase, const int &indexBase, const int &nElements, const int &nNodes, const double *nodeArray, const int *elementNodesArray, const int *nodeMaterialTypes, const int *elementMaterialTypes)
bool writeTetgenMeshNodesAndElements(const char *filebase, const int &indexBase, const int &nElements, const int &nNodes, const double *nodeArray, const int *elementNodesArray, const int *nodeMaterialTypes, const int *elementMaterialTypes)
bool write2dmMeshNodesAndElements(const char *filebase, const int &indexBase, const int &nElements, const int &nNodes, const double *nodeArray, const int *elementNodesArray, const int *elementMaterialTypes)
bool readTetgenMeshNodesAndElements(const char *filebase, const int &indexBase, int &nElements, int &nNodes, std::vector< double > &nodeArray, std::vector< int > &elementNodesArray, std::vector< int > &nodeMaterialTypes, std::vector< int > &elementMaterialTypes, const int &defaultElementMaterialType, const int &defaultNodeMaterialType)
bool readTriangleElementBoundaries(const char *filebase, const int &indexBase, bool &hasMarkers, int &nElementBoundaries, std::vector< int > &elementBoundaryNodesArray, std::vector< int > &elementBoundaryMaterialTypesArray, const int &defaultBoundaryMaterialType)
bool writeTetgenElementBoundaryNodes(const char *filebase, const int &indexBase, const int &nElementBoundariesToWrite, const int *elementBoundaryNodesArray, const int *elementBoundaryMaterialTypes, const bool &writeExteriorElementBoundariesOnly, const int *exteriorElementBoundariesArray)
int * elementBoundaryNodesArray
double * elementBoundaryBarycentersArray
int nNodes_elementBoundary
double * elementDiametersArray
int * elementBoundaryElementsArray
double * nodeSupportArray
double * elementBoundaryDiametersArray
int * elementNeighborsArray
int * exteriorElementBoundariesArray
int * elementBoundaryLocalElementBoundariesArray
int nExteriorElementBoundaries_global
int nElementBoundaries_element
int * elementBoundaryMaterialTypes
int max_nNodeNeighbors_node
int * elementMaterialTypes
double * elementBarycentersArray
double * elementInnerDiametersArray
int nInteriorElementBoundaries_global
int * elementBoundariesArray
int * interiorElementBoundariesArray
double * nodeDiametersArray
int nElementBoundaries_global
int ** elementParentsArray
int ** elementChildrenArray
int ** elementChildrenOffsets