9 #pragma mark * local ( static ) function prototypes *
11 static double CurrentTime(
void)
26 #pragma mark * exported function implementations *
33 #include PROTEUS_TRIANGLE_H
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 register 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);
633 for(
int ebN=0;ebN<2;ebN++)
635 register int nodes[1];
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++)
793 register int nodes[2];
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++)
962 register int nodes[2];
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();
1122 using namespace std;
1126 ElementBoundaryElementsType elementBoundaryElements;
1127 start=CurrentTime();
1132 register int nodes[3];
1137 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1139 elementBoundaryElements[ebt].right=eN;
1140 elementBoundaryElements[ebt].right_ebN_element=ebN;
1144 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1147 stop = CurrentTime();
1153 start = CurrentTime();
1154 typedef set<int> ElementBoundariesType;
1155 ElementBoundariesType interiorElementBoundaries,exteriorElementBoundaries;
1161 stop = CurrentTime();
1164 start = CurrentTime();
1167 eb != elementBoundaryElements.end();
1179 if(eb->second.right != -1)
1181 interiorElementBoundaries.insert(ebN);
1185 exteriorElementBoundaries.insert(ebN);
1187 if (eb->second.right != -1)
1197 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1199 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1206 typedef set<NodeTuple<2> > EdgesType;
1222 EdgesType::iterator edge_p=edges.begin();
1223 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
1232 typedef vector<set<int> > NodeStarType;
1245 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1250 stop = CurrentTime();
1254 typedef vector<set<int> > NodeElementsStarType;
1268 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1310 using namespace std;
1314 start=CurrentTime();
1316 int lface[6][4] = {{0,1,2,3},
1327 register int nodes[4];
1334 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1336 elementBoundaryElements[ebt].right=eN;
1337 elementBoundaryElements[ebt].right_ebN_element=ebN;
1341 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1344 stop = CurrentTime();
1350 start = CurrentTime();
1351 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1358 stop = CurrentTime();
1362 start = CurrentTime();
1365 eb != elementBoundaryElements.end();
1378 if(eb->second.right != -1)
1380 interiorElementBoundaries.insert(ebN);
1384 exteriorElementBoundaries.insert(ebN);
1387 if (eb->second.right != -1)
1397 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1399 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1401 set<NodeTuple<2> > edges;
1403 int ledge[12][2] = {{0,1},{1,2},{2,3},{3,0},
1404 {0,4},{1,5},{2,6},{3,7},
1405 {4,5},{5,6},{6,7},{7,4}};
1414 for (
int e=0;e<12;e++)
1427 set<NodeTuple<2> >::iterator edge_p=edges.begin();
1428 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
1451 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1453 stop = CurrentTime();
1473 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1516 using namespace std;
1520 int n2 = (mesh.
px+1)*(mesh.
py+1)-1 ;
1521 int n3 = (mesh.
px+1)*mesh.
py;
1524 int n4 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n0;
1525 int n5 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n1;
1526 int n6 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n2;
1527 int n7 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n3;
1533 using namespace std;
1537 start=CurrentTime();
1539 int lface[6][4] = {{n0,n1,n2,n3},
1552 register int nodes[4];
1559 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1561 elementBoundaryElements[ebt].right=eN;
1562 elementBoundaryElements[ebt].right_ebN_element=ebN;
1566 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1569 stop = CurrentTime();
1575 start = CurrentTime();
1576 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1583 stop = CurrentTime();
1587 start = CurrentTime();
1590 eb != elementBoundaryElements.end();
1603 if(eb->second.right != -1)
1605 interiorElementBoundaries.insert(ebN);
1609 exteriorElementBoundaries.insert(ebN);
1612 if (eb->second.right != -1)
1622 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1624 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1626 set<NodeTuple<2> > edges;
1631 int ledge[12][2] = {{n0,n1},{n1,n2},{n2,n3},{n3,n0},
1632 {n0,n4},{n1,n5},{n2,n6},{n3,n7},
1633 {n4,n5},{n5,n6},{n6,n7},{n7,n4}};
1642 for (
int e=0;e<12;e++)
1655 set<NodeTuple<2> >::iterator edge_p=edges.begin();
1656 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
1676 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1678 stop = CurrentTime();
1697 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1741 int n2 = (mesh.
px+1)*(mesh.
py+1)-1 ;
1742 int n3 = (mesh.
px+1)*mesh.
py;
1745 int n4 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n0;
1746 int n5 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n1;
1747 int n6 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n2;
1748 int n7 = (mesh.
px+1)*(mesh.
py+1)*mesh.
pz + n3;
1753 using namespace std;
1758 int> elementBoundaryIds;
1759 start=CurrentTime();
1761 int lface[6][4] = {{n0,n1,n2,n3},
1773 register int nodes[4];
1779 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1781 elementBoundaryElements[ebt].right=eN;
1782 elementBoundaryElements[ebt].right_ebN_element=ebN;
1788 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1789 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
1792 stop = CurrentTime();
1798 start = CurrentTime();
1799 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1804 stop = CurrentTime();
1808 start = CurrentTime();
1810 eb != elementBoundaryElements.end();
1813 int ebN = elementBoundaryIds[eb->first];
1824 if(eb->second.right != -1)
1826 interiorElementBoundaries.insert(ebN);
1830 exteriorElementBoundaries.insert(ebN);
1832 if (eb->second.right != -1)
1842 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
1844 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
1884 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
1886 stop = CurrentTime();
1904 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
1947 using namespace std;
1953 int> elementBoundaryIds;
1954 start=CurrentTime();
1956 for(
int ebN=0;ebN<2;ebN++)
1959 register int nodes[1];
1962 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
1964 elementBoundaryElements[ebt].right=eN;
1965 elementBoundaryElements[ebt].right_ebN_element=ebN;
1966 assert(elementBoundaryIds[ebt] == ebN_global);
1970 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
1971 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
1974 stop = CurrentTime();
1980 start = CurrentTime();
1981 set<int> interiorElementBoundaries,exteriorElementBoundaries;
1986 stop = CurrentTime();
1990 start = CurrentTime();
1992 eb != elementBoundaryElements.end();
1995 int ebN = elementBoundaryIds[eb->first];
2002 if(eb->second.right != -1)
2004 interiorElementBoundaries.insert(ebN);
2008 exteriorElementBoundaries.insert(ebN);
2010 if (eb->second.right != -1)
2020 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2022 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2044 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2049 stop = CurrentTime();
2064 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2105 using namespace std;
2114 int> elementBoundaryIds;
2115 start=CurrentTime();
2117 for(
int ebN=0;ebN<3;ebN++)
2120 register int nodes[2];
2124 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2126 elementBoundaryElements[ebt].right=eN;
2127 elementBoundaryElements[ebt].right_ebN_element=ebN;
2128 assert(elementBoundaryIds[ebt] == ebN_global);
2132 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2133 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2136 stop = CurrentTime();
2142 start = CurrentTime();
2143 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2148 stop = CurrentTime();
2152 start = CurrentTime();
2154 eb != elementBoundaryElements.end();
2157 int ebN = elementBoundaryIds[eb->first];
2165 if(eb->second.right != -1)
2167 interiorElementBoundaries.insert(ebN);
2171 exteriorElementBoundaries.insert(ebN);
2173 if (eb->second.right != -1)
2183 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2185 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2187 set<NodeTuple<2> > edges;
2202 for (set<
NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
2219 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2221 stop = CurrentTime();
2239 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2279 using namespace std;
2288 int> elementBoundaryIds;
2289 start=CurrentTime();
2291 for(
int ebN=0;ebN<4;ebN++)
2294 register int nodes[2];
2298 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2300 elementBoundaryElements[ebt].right=eN;
2301 elementBoundaryElements[ebt].right_ebN_element=ebN;
2302 assert(elementBoundaryIds[ebt] == ebN_global);
2306 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2307 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2310 stop = CurrentTime();
2316 start = CurrentTime();
2317 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2322 stop = CurrentTime();
2326 start = CurrentTime();
2328 eb != elementBoundaryElements.end();
2331 int ebN = elementBoundaryIds[eb->first];
2339 if(eb->second.right != -1)
2341 interiorElementBoundaries.insert(ebN);
2345 exteriorElementBoundaries.insert(ebN);
2347 if (eb->second.right != -1)
2357 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2359 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2361 set<NodeTuple<2> > edges;
2375 for (set<
NodeTuple<2> >::iterator edgeTuple_p=edges.begin();edgeTuple_p != edges.end();edgeTuple_p++,edgeN++)
2392 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2394 stop = CurrentTime();
2412 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2455 using namespace std;
2460 int> elementBoundaryIds;
2461 start=CurrentTime();
2467 register int nodes[3];
2472 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2474 elementBoundaryElements[ebt].right=eN;
2475 elementBoundaryElements[ebt].right_ebN_element=ebN;
2476 assert(elementBoundaryIds[ebt] == ebN_global);
2480 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2481 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2484 stop = CurrentTime();
2490 start = CurrentTime();
2491 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2496 stop = CurrentTime();
2500 start = CurrentTime();
2502 eb != elementBoundaryElements.end();
2505 int ebN = elementBoundaryIds[eb->first];
2515 if(eb->second.right != -1)
2517 interiorElementBoundaries.insert(ebN);
2521 exteriorElementBoundaries.insert(ebN);
2523 if (eb->second.right != -1)
2533 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2535 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2537 set<NodeTuple<2> > edges;
2552 set<NodeTuple<2> >::iterator edge_p=edges.begin();
2553 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
2571 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2573 stop = CurrentTime();
2591 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2635 using namespace std;
2641 int> elementBoundaryIds;
2642 start=CurrentTime();
2644 for(
int ebN=0;ebN<2;ebN++)
2647 register int nodes[1];
2650 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2652 elementBoundaryElements[ebt].right=eN;
2653 elementBoundaryElements[ebt].right_ebN_element=ebN;
2654 assert(elementBoundaryIds[ebt] == ebN_global);
2658 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2659 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2662 stop = CurrentTime();
2668 start = CurrentTime();
2669 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2674 stop = CurrentTime();
2678 start = CurrentTime();
2680 eb != elementBoundaryElements.end();
2683 int ebN = elementBoundaryIds[eb->first];
2690 if(eb->second.right != -1)
2692 interiorElementBoundaries.insert(ebN);
2696 exteriorElementBoundaries.insert(ebN);
2698 if (eb->second.right != -1)
2708 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2710 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2726 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2731 stop = CurrentTime();
2746 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2787 using namespace std;
2796 int> elementBoundaryIds;
2797 start=CurrentTime();
2799 for(
int ebN=0;ebN<3;ebN++)
2802 register int nodes[2];
2806 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2808 elementBoundaryElements[ebt].right=eN;
2809 elementBoundaryElements[ebt].right_ebN_element=ebN;
2810 assert(elementBoundaryIds[ebt] == ebN_global);
2814 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2815 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2818 stop = CurrentTime();
2824 start = CurrentTime();
2825 set<int> interiorElementBoundaries,exteriorElementBoundaries;
2830 stop = CurrentTime();
2834 start = CurrentTime();
2836 eb != elementBoundaryElements.end();
2839 int ebN = elementBoundaryIds[eb->first];
2847 if(eb->second.right != -1)
2849 interiorElementBoundaries.insert(ebN);
2853 exteriorElementBoundaries.insert(ebN);
2855 if (eb->second.right != -1)
2865 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
2867 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
2902 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
2904 stop = CurrentTime();
2922 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
2962 using namespace std;
2971 int> elementBoundaryIds;
2972 start=CurrentTime();
2974 for(
int ebN=0;ebN<4;ebN++)
2977 register int nodes[2];
2981 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
2983 elementBoundaryElements[ebt].right=eN;
2984 elementBoundaryElements[ebt].right_ebN_element=ebN;
2985 assert(elementBoundaryIds[ebt] == ebN_global);
2989 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
2990 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
2993 stop = CurrentTime();
2999 start = CurrentTime();
3000 set<int> interiorElementBoundaries,exteriorElementBoundaries;
3005 stop = CurrentTime();
3009 start = CurrentTime();
3011 eb != elementBoundaryElements.end();
3014 int ebN = elementBoundaryIds[eb->first];
3022 if(eb->second.right != -1)
3024 interiorElementBoundaries.insert(ebN);
3028 exteriorElementBoundaries.insert(ebN);
3030 if (eb->second.right != -1)
3040 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
3042 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
3057 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3059 stop = CurrentTime();
3077 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3120 using namespace std;
3125 int> elementBoundaryIds;
3126 start=CurrentTime();
3132 register int nodes[3];
3138 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
3140 elementBoundaryElements[ebt].right=eN;
3141 elementBoundaryElements[ebt].right_ebN_element=ebN;
3142 assert(elementBoundaryIds[ebt] == ebN_global);
3146 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
3147 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
3150 stop = CurrentTime();
3156 start = CurrentTime();
3157 set<int> interiorElementBoundaries,exteriorElementBoundaries;
3162 stop = CurrentTime();
3166 start = CurrentTime();
3168 eb != elementBoundaryElements.end();
3171 int ebN = elementBoundaryIds[eb->first];
3181 if(eb->second.right != -1)
3183 interiorElementBoundaries.insert(ebN);
3187 exteriorElementBoundaries.insert(ebN);
3189 if (eb->second.right != -1)
3199 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
3201 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
3238 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3257 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3299 using namespace std;
3304 int> elementBoundaryIds;
3305 start=CurrentTime();
3307 int lface[6][4] = {{0,1,2,3},
3319 register int nodes[4];
3325 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
3327 elementBoundaryElements[ebt].right=eN;
3328 elementBoundaryElements[ebt].right_ebN_element=ebN;
3332 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
3333 elementBoundaryIds.insert(elementBoundaryIds.end(),make_pair(ebt,ebN_global));
3336 stop = CurrentTime();
3342 start = CurrentTime();
3343 set<int> interiorElementBoundaries,exteriorElementBoundaries;
3348 stop = CurrentTime();
3352 start = CurrentTime();
3354 eb != elementBoundaryElements.end();
3357 int ebN = elementBoundaryIds[eb->first];
3368 if(eb->second.right != -1)
3370 interiorElementBoundaries.insert(ebN);
3374 exteriorElementBoundaries.insert(ebN);
3376 if (eb->second.right != -1)
3386 for (set<int>::iterator ebN=interiorElementBoundaries.begin();ebN != interiorElementBoundaries.end(); ebN++,ebNI++)
3388 for (set<int>::iterator ebN=exteriorElementBoundaries.begin();ebN != exteriorElementBoundaries.end(); ebN++,ebNE++)
3426 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
3428 stop = CurrentTime();
3446 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
3486 register double dx,dy,dz;
3487 dx = nodeArray[nL*3+0] - nodeArray[nR*3+0];
3488 dy = nodeArray[nL*3+1] - nodeArray[nR*3+1];
3489 dz = nodeArray[nL*3+2] - nodeArray[nR*3+2];
3490 return sqrt(dx*dx+dy*dy+dz*dz);
3493 inline double triangleArea(
int n0,
int n1,
int n2,
const double* nodeArray)
3495 register double va[3],vb[3];
3496 va[0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
3497 va[1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
3498 va[2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
3499 vb[0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
3500 vb[1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
3501 vb[2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
3502 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]));
3507 register double t[3][3];
3508 t[0][0] = nodeArray[n1*3+0] - nodeArray[n0*3+0];
3509 t[0][1] = nodeArray[n1*3+1] - nodeArray[n0*3+1];
3510 t[0][2] = nodeArray[n1*3+2] - nodeArray[n0*3+2];
3512 t[1][0] = nodeArray[n2*3+0] - nodeArray[n0*3+0];
3513 t[1][1] = nodeArray[n2*3+1] - nodeArray[n0*3+1];
3514 t[1][2] = nodeArray[n2*3+2] - nodeArray[n0*3+2];
3516 t[2][0] = nodeArray[n3*3+0] - nodeArray[n0*3+0];
3517 t[2][1] = nodeArray[n3*3+1] - nodeArray[n0*3+1];
3518 t[2][2] = nodeArray[n3*3+2] - nodeArray[n0*3+2];
3519 return fabs(t[0][0]*(t[1][1]*t[2][2] - t[1][2]*t[2][1]) - \
3520 t[0][1]*(t[1][0]*t[2][2] - t[1][2]*t[2][0]) + \
3521 t[0][2]*(t[1][0]*t[2][1] - t[1][1]*t[2][0]))/6.0;
3575 double volume, surfaceArea=0.0, hMin=0.0,hMax=0.0;
3582 for (
int ebN=0;ebN<4;ebN++)
3589 hMin = 6.0*volume/surfaceArea;
3601 mesh.
h = fmax(hMax,mesh.
h);
3629 inline double hexahedronVolume(
int n0,
int n1,
int n2,
int n3,
int n4,
int n5,
int n6,
int n7,
const double* nodeArray)
3631 register double t[3];
3632 t[0] = nodeArray[n0*3+0] - nodeArray[n1*3+0];
3633 t[1] = nodeArray[n0*3+1] - nodeArray[n3*3+1];
3634 t[2] = nodeArray[n0*3+2] - nodeArray[n4*3+2];
3636 return fabs(t[0]*t[1]*t[2]);
3679 double volume, hMin=0.0,hMax=0.0;
3709 mesh.
h = fmax(hMax,mesh.
h);
3776 double area=0.0, perimeter=0.0,h=0.0, hMin=0.0,hMax=0.0;
3789 hMax = fmax(hMax,h);
3791 hMin = 4.0*area/perimeter;
3796 mesh.
h = fmax(hMax,mesh.
h);
3872 double area=0.0, perimeter=0.0,h=0.0, hMin=1.0e6,hMax=0.0;
3890 hMax = fmax(hMax,h);
3891 hMin = fmin(hMin,h);
3897 mesh.
h = fmax(hMax,mesh.
h);
4016 using namespace std;
4017 multilevelMesh.
nLevels = nLevels;
4019 for(
int i=0;i<nLevels;i++)
4025 for(
int i=1;i<nLevels;i++)
4028 set<Node> newNodeSet;
4029 set<Node>::iterator nodeItr;
4030 pair<set<Node>::iterator,
bool> ret;
4054 midpoints[0].
nN = nN_new;
4058 else if (averageNewNodeFlags)
4064 newNodeSet.insert(midpoints[0]);
4084 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4101 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4111 using namespace std;
4112 multilevelMesh.
nLevels = nLevels;
4114 for(
int i=0;i<nLevels;i++)
4120 for(
int i=1;i<nLevels;i++)
4123 set<Node> newNodeSet;
4124 set<Node>::iterator nodeItr;
4125 pair<set<Node>::iterator,
bool> ret;
4152 for(
int nN_element_0=0,nN_midpoint=0;nN_element_0<mesh.
nNodes_element;nN_element_0++)
4153 for(
int nN_element_1=nN_element_0+1;nN_element_1<mesh.
nNodes_element;nN_element_1++,nN_midpoint++)
4157 midpoints[nN_midpoint]);
4158 nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
4159 if(nodeItr == newNodeSet.end())
4161 midpoints[nN_midpoint].
nN = nN_new;
4165 else if (averageNewNodeFlags)
4171 newNodeSet.insert(midpoints[nN_midpoint]);
4175 midpoints[nN_midpoint].
nN = nodeItr->nN;
4205 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4222 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4236 using namespace std;
4237 multilevelMesh.
nLevels = nLevels;
4239 for(
int i=0;i<nLevels;i++)
4245 for(
int i=1;i<nLevels;i++)
4247 std::cout<<
"Quad refinement not imlemented "<<i<<endl;
4255 using namespace std;
4256 multilevelMesh.
nLevels = nLevels;
4258 for(
int i=0;i<nLevels;i++)
4264 for(
int i=1;i<nLevels;i++)
4267 set<Node> newNodeSet;
4268 set<Node>::iterator nodeItr;
4269 pair<set<Node>::iterator,
bool> ret;
4311 for(
int nN_element_0=0,nN_midpoint=0;nN_element_0<mesh.
nNodes_element;nN_element_0++)
4312 for(
int nN_element_1=nN_element_0+1;nN_element_1<mesh.
nNodes_element;nN_element_1++,nN_midpoint++)
4316 midpoints[nN_midpoint]);
4317 nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
4318 if(nodeItr == newNodeSet.end())
4320 midpoints[nN_midpoint].
nN = nN_new;
4324 else if (averageNewNodeFlags)
4330 newNodeSet.insert(midpoints[nN_midpoint]);
4334 midpoints[nN_midpoint].
nN = nodeItr->nN;
4360 for(
int dN=1;dN<3;dN++)
4362 dd =
edgeLength(midpoints[dN],midpoints[5-dN]);
4368 else if (dd == mind)
4370 if(midpoints[dN] < midpoints[mindN])
4372 else if (!(midpoints[mindN] < midpoints[dN]) && (midpoints[5-dN] < midpoints[5-mindN]))
4399 else if (mindN == 1)
4456 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4530 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
4540 using namespace std;
4541 multilevelMesh.
nLevels = nLevels;
4543 for(
int i=0;i<nLevels;i++)
4549 for(
int i=1;i<nLevels;i++)
4551 std::cout<<
"Hexahedron refinement not implemented"<<std::endl;
4558 const int& nSpace_global)
4566 int eN_left_parent = levelElementParentsArray[eN_left];
4567 int eN_right_parent= levelElementParentsArray[eN_right];
4568 if (eN_left_parent == eN_right_parent)
4573 int left_parent_ebN = -1;
4575 ebN_left_parent_element++)
4578 ebN_left_parent_element]
4581 left_parent_ebN = ebN_left_parent_element;
4595 int eN_parent = levelElementParentsArray[eN];
4596 int lastExteriorElementBoundaryOnParent=-1;
4597 int nExteriorElementBoundariesOnParent =0;
4602 lastExteriorElementBoundaryOnParent = ebN_element;
4603 nExteriorElementBoundariesOnParent++;
4606 assert(nExteriorElementBoundariesOnParent > 0); assert(lastExteriorElementBoundaryOnParent >= 0);
4607 if (nExteriorElementBoundariesOnParent > 1)
4611 if (nSpace_global == 1)
4618 const double n_child = (childMesh.
nodeArray[nN_child1*3+0]-childMesh.
nodeArray[nN_child0*3+0])/
4623 const double n_parent0 = (parentMesh.
nodeArray[nN_parent1*3+0]-parentMesh.
nodeArray[nN_parent0*3+0])/
4625 const double n_parent1 = (parentMesh.
nodeArray[nN_parent0*3+0]-parentMesh.
nodeArray[nN_parent1*3+0])/
4628 int ebN_parent = -1;
4629 if (fabs(n_child-n_parent1) < 1.0e-8)
4631 else if (fabs(n_child-n_parent0) < 1.0e-8)
4633 assert(ebN_parent >= 0);
4636 else if (nSpace_global == 2)
4640 double n_child[2] = {childMesh.
nodeArray[nN_ebN_child[1]*3+1]-childMesh.
nodeArray[nN_ebN_child[0]*3+1],
4642 double tmp = sqrt(n_child[0]*n_child[0] + n_child[1]*n_child[1]);
4644 n_child[0] /= tmp; n_child[1] /= tmp;
4646 int ebN_parent = -1;
4653 double n_parent[2] = {parentMesh.
nodeArray[nN_ebN_parent[1]*3+1]-parentMesh.
nodeArray[nN_ebN_parent[0]*3+1],
4654 parentMesh.
nodeArray[nN_ebN_parent[0]*3+0]-parentMesh.
nodeArray[nN_ebN_parent[1]*3+0]};
4655 tmp = sqrt(n_parent[0]*n_parent[0] + n_parent[1]*n_parent[1]);
4657 n_parent[0] /= tmp; n_parent[1] /= tmp;
4660 tmp = n_parent[0]*n_child[0] + n_parent[1]*n_child[1];
4662 if (fabs(sqrt(tmp*tmp) - 1.0) < 1.0e-8)
4664 ebN_parent = ebN_cur;
4668 assert(ebN_parent >= 0);
4676 double n_child[3] = {0.,0.,0.};
4678 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])
4679 -(nN_ebN_child_1_x[2]-nN_ebN_child_0_x[2])*(nN_ebN_child_2_x[1]-nN_ebN_child_0_x[1]);
4681 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])
4682 -(nN_ebN_child_1_x[0]-nN_ebN_child_0_x[0])*(nN_ebN_child_2_x[2]-nN_ebN_child_0_x[2]);
4684 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])
4685 -(nN_ebN_child_1_x[1]-nN_ebN_child_0_x[1])*(nN_ebN_child_2_x[0]-nN_ebN_child_0_x[0]);
4687 double tmp = sqrt(n_child[0]*n_child[0] + n_child[1]*n_child[1] + n_child[2]*n_child[2]);
4689 n_child[0] /= tmp; n_child[1] /= tmp; n_child[2] /= tmp;
4691 int ebN_parent = -1;
4692 double n_parent[3] = {0.,0.,0.};
4701 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])
4702 -(nN_ebN_parent_1_x[2]-nN_ebN_parent_0_x[2])*(nN_ebN_parent_2_x[1]-nN_ebN_parent_0_x[1]);
4704 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])
4705 -(nN_ebN_parent_1_x[0]-nN_ebN_parent_0_x[0])*(nN_ebN_parent_2_x[2]-nN_ebN_parent_0_x[2]);
4707 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])
4708 -(nN_ebN_parent_1_x[1]-nN_ebN_parent_0_x[1])*(nN_ebN_parent_2_x[0]-nN_ebN_parent_0_x[0]);
4710 tmp = sqrt(n_parent[0]*n_parent[0] + n_parent[1]*n_parent[1] + n_parent[2]*n_parent[2]);
4712 n_parent[0] /= tmp; n_parent[1] /= tmp; n_parent[2] /= tmp;
4716 tmp = n_parent[0]*n_child[0] + n_parent[1]*n_child[1] + n_parent[2]*n_child[2];
4718 if (fabs(sqrt(tmp*tmp) - 1.0) < 1.0e-8)
4720 ebN_parent = ebN_cur;
4725 assert(ebN_parent >= 0);
4733 lastExteriorElementBoundaryOnParent];
4745 using namespace std;
4748 string word,elementType;
4753 if(word ==
"MESH1D")
4755 elementType =
"E2E";
4759 else if(word ==
"MESH2D")
4761 elementType =
"E3T";
4765 else if (word ==
"MESH3D")
4767 elementType =
"E4T";
4773 cerr<<
"Unrecognized mesh type"<<endl;
4776 vector<vector<int> > elementNodesVector;
4778 vector<int> materialTypes;
4781 while(word == elementType)
4786 meshFile>>nodes[nN];
4789 elementNodesVector.push_back(nodes);
4791 materialTypes.push_back(type);
4809 using namespace std;
4811 meshFile<<
"ND"<<setw(7)<<nN+1
4812 <<scientific<<setprecision(8)<<setw(16)<<mesh.
nodeArray[3*nN+0]
4813 <<scientific<<setprecision(8)<<setw(16)<<mesh.
nodeArray[3*nN+1]
4814 <<scientific<<setprecision(8)<<setw(16)<<mesh.
nodeArray[3*nN+2]
4822 meshFile<<
"Try to write something"<<std::endl;
4823 using namespace std;
4829 elementType =
"E2E";
4830 meshFile<<
"MESH1D"<<endl;
4835 elementType =
"E3T";
4836 meshFile<<
"MESH2D"<<endl;
4841 elementType =
"E4T";
4842 meshFile<<
"MESH3D"<<endl;
4846 cerr<<
"Unknown element type"<<endl;
4851 meshFile<<elementType;
4852 meshFile<<setw(width)<<eN+1;
4866 assert(trimesh); assert(trimesh->trianglelist);
4875 if (trimesh->triangleattributelist)
4877 for (
int eN = 0; eN < trimesh->numberoftriangles; eN++)
4878 mesh.
elementMaterialTypes[eN] = trimesh->triangleattributelist[eN*trimesh->numberoftriangleattributes+0];
4889 trimesh->trianglelist[eN*trimesh->numberofcorners+ebN]-base;
4900 assert(trimesh); assert(trimesh->pointlist);
4913 if (trimesh->pointmarkerlist)
4915 for (
int nN = 0; nN < trimesh->numberofpoints; nN++)
4925 trimesh->pointlist[nN*2+0];
4927 trimesh->pointlist[nN*2+1];
4939 if (!trimesh->edgelist)
4941 if (!trimesh->edgemarkerlist)
4948 std::map<NodeTuple<2>,
int> triangleEdgeMarker;
4949 for (
int ebN = 0; ebN < trimesh->numberofedges; ebN++)
4952 nodes[0] = trimesh->edgelist[ebN*2+0]-base;
4953 nodes[1] = trimesh->edgelist[ebN*2+1]-base;
4955 triangleEdgeMarker[ebt] = trimesh->edgemarkerlist[ebN];
4978 using namespace std;
4982 int nElements_file, nNodes_file;
4983 vector<double> nodeArray_file;
4984 vector<int> elementNodesArray_file;
4985 vector<int> nodeMaterialTypes_file;
4986 vector<int> elementMaterialTypes_file;
4992 elementNodesArray_file,
4993 nodeMaterialTypes_file,
4994 elementMaterialTypes_file,
5010 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5013 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5022 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5024 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5036 using namespace std;
5070 using namespace std;
5074 int nElementBoundaries_file;
5075 vector<int> elementBoundaryNodesArray_file;
5076 vector<int> elementBoundaryMaterialTypes_file;
5077 bool elementBoundaryMaterialTypesInFile =
false;
5080 elementBoundaryMaterialTypesInFile,
5081 nElementBoundaries_file,
5082 elementBoundaryNodesArray_file,
5083 elementBoundaryMaterialTypes_file,
5090 if (!elementBoundaryMaterialTypesInFile)
5096 map<NodeTuple<2>,
int> triangleElementBoundaryMaterialTypes;
5098 for (
int ebN = 0; ebN < nElementBoundaries_file; ebN++)
5100 register int nodes[2];
5101 nodes[0] = elementBoundaryNodesArray_file[ebN*2+0];
5102 nodes[1] = elementBoundaryNodesArray_file[ebN*2+1];
5104 triangleElementBoundaryMaterialTypes[ttuple] = elementBoundaryMaterialTypes_file[ebN];
5108 register int nodes[2];
5126 using namespace std;
5130 int nElements_file, nNodes_file;
5131 vector<double> nodeArray_file;
5132 vector<int> elementNodesArray_file;
5133 vector<int> nodeMaterialTypes_file;
5134 vector<int> elementMaterialTypes_file;
5140 elementNodesArray_file,
5141 nodeMaterialTypes_file,
5142 elementMaterialTypes_file,
5157 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5161 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5170 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5172 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5187 using namespace std;
5191 int nElementBoundaries_file;
5192 vector<int> elementBoundaryNodesArray_file;
5193 vector<int> elementBoundaryMaterialTypes_file;
5194 bool elementBoundaryMaterialTypesInFile =
false;
5197 elementBoundaryMaterialTypesInFile,
5198 nElementBoundaries_file,
5199 elementBoundaryNodesArray_file,
5200 elementBoundaryMaterialTypes_file,
5207 if (!elementBoundaryMaterialTypesInFile)
5213 using namespace std;
5218 start=CurrentTime();
5223 register int nodes[3];
5228 if(elementBoundaryElements.find(ebt) != elementBoundaryElements.end())
5230 elementBoundaryElements[ebt].right=eN;
5231 elementBoundaryElements[ebt].right_ebN_element=ebN;
5235 elementBoundaryElements.insert(elementBoundaryElements.end(),make_pair(ebt,
ElementNeighbors(eN,ebN)));
5238 stop = CurrentTime();
5244 start = CurrentTime();
5252 stop = CurrentTime();
5256 start = CurrentTime();
5259 eb != elementBoundaryElements.end();
5280 if (eb->second.right != -1)
5296 set<NodeTuple<2> > edges;
5311 set<NodeTuple<2> >::iterator edge_p=edges.begin();
5312 for (
int edgeN=0;edgeN<int(edges.size());edgeN++)
5331 for (set<int>::iterator nN_star=nodeStar[nN].begin();nN_star!=nodeStar[nN].end();nN_star++,offset++)
5333 stop = CurrentTime();
5352 for (set<int>::iterator eN_star = nodeElementsStar[nN].begin(); eN_star != nodeElementsStar[nN].end();
5400 map<NodeTuple<3>,
int> tetgenElementBoundaryMaterialTypes;
5402 for (
int ebN = 0; ebN < nElementBoundaries_file; ebN++)
5404 register int nodes[3];
5405 nodes[0] = elementBoundaryNodesArray_file[ebN*3+0];
5406 nodes[1] = elementBoundaryNodesArray_file[ebN*3+1];
5407 nodes[2] = elementBoundaryNodesArray_file[ebN*3+2];
5409 tetgenElementBoundaryMaterialTypes[ttuple] = elementBoundaryMaterialTypes_file[ebN];
5413 register int nodes[3];
5426 map<NodeTuple<3>,
int> tetgenElementBoundaryMaterialTypes;
5428 for (
int ebNE = 0; ebNE < nElementBoundaries_file; ebNE++)
5430 register int nodes[3];
5431 nodes[0] = elementBoundaryNodesArray_file[ebNE*3+0];
5432 nodes[1] = elementBoundaryNodesArray_file[ebNE*3+1];
5433 nodes[2] = elementBoundaryNodesArray_file[ebNE*3+2];
5435 tetgenElementBoundaryMaterialTypes[ttuple] = elementBoundaryMaterialTypes_file[ebNE];
5440 register int nodes[3];
5461 using namespace std;
5478 bool writeExteriorElementBoundariesOnly =
true;
5482 nElementBoundariesToWrite,
5485 writeExteriorElementBoundariesOnly,
5498 using namespace std;
5501 std::string meshFilename= std::string(filebase)+
".3dm";
5502 std::ifstream meshFile(meshFilename.c_str());
5503 if (!meshFile.good())
5505 std::cerr<<
"read3DM cannot open file "
5506 <<meshFilename<<std::endl;
5511 std::string fileType;
5513 if (fileType !=
"MESH3D")
5515 std::cerr<<
"read3DM does not recognize filetype "
5516 <<fileType<<std::endl;
5520 std::string firstWord;
5521 meshFile>>firstWord;
5523 std::vector<int> elementNodesArray_file;
5524 std::vector<int> elementMaterialTypes_file;
5525 int eN,n0,n1,n2,n3,emt;
5526 while (firstWord ==
"E4T")
5529 meshFile>>eN>>n0>>n1>>n2>>n3>>emt;
5530 elementNodesArray_file.push_back(n0-indexBase);
5531 elementNodesArray_file.push_back(n1-indexBase);
5532 elementNodesArray_file.push_back(n2-indexBase);
5533 elementNodesArray_file.push_back(n3-indexBase);
5534 elementMaterialTypes_file.push_back(emt-indexBase);
5535 meshFile>>firstWord;
5537 std::vector<int> nodeMaterialTypes_file;
5538 std::vector<double> nodeArray_file;
5542 while (!meshFile.eof() && firstWord ==
"ND")
5545 meshFile>>nN>>x>>y>>
z;
5546 nodeArray_file.push_back(x);
5547 nodeArray_file.push_back(y);
5548 nodeArray_file.push_back(
z);
5549 nodeMaterialTypes_file.push_back(0);
5550 meshFile>>firstWord;
5555 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5559 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5564 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5569 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5582 using namespace std;
5585 std::string meshFilename= std::string(filebase)+
".3dm";
5586 std::ifstream meshFile(meshFilename.c_str());
5587 if (!meshFile.good())
5589 std::cerr<<
"read2DM cannot open file "
5590 <<meshFilename<<std::endl;
5595 std::string fileType;
5597 if (fileType !=
"MESH2D")
5599 std::cerr<<
"read2DM does not recognize filetype "
5600 <<fileType<<std::endl;
5604 std::string firstWord;
5605 meshFile>>firstWord;
5607 std::vector<int> elementNodesArray_file;
5608 std::vector<int> elementMaterialTypes_file;
5609 int eN,n0,n1,n2,emt;
5610 while (firstWord ==
"E3T")
5613 meshFile>>eN>>n0>>n1>>n2>>emt;
5614 elementNodesArray_file.push_back(n0-indexBase);
5615 elementNodesArray_file.push_back(n1-indexBase);
5616 elementNodesArray_file.push_back(n2-indexBase);
5617 elementMaterialTypes_file.push_back(emt-indexBase);
5618 meshFile>>firstWord;
5620 std::vector<int> nodeMaterialTypes_file;
5621 std::vector<double> nodeArray_file;
5625 while (!meshFile.eof() && firstWord ==
"ND")
5628 meshFile>>nN>>x>>y>>
z;
5629 nodeArray_file.push_back(x);
5630 nodeArray_file.push_back(y);
5631 nodeArray_file.push_back(
z);
5632 nodeMaterialTypes_file.push_back(0);
5633 meshFile>>firstWord;
5638 copy(nodeArray_file.begin(),nodeArray_file.end(),mesh.
nodeArray);
5642 copy(nodeMaterialTypes_file.begin(),nodeMaterialTypes_file.end(),
5647 copy(elementNodesArray_file.begin(),elementNodesArray_file.end(),
5652 copy(elementMaterialTypes_file.begin(),elementMaterialTypes_file.end(),
5665 using namespace std;
5668 std::string meshFilename= std::string(filebase)+
".mesh";
5669 std::ifstream meshFile(meshFilename.c_str());
5673 if (!meshFile.good())
5675 std::cerr<<
"readHex cannot open file "
5676 <<meshFilename<<std::endl;
5681 std::string fileType;
5683 if (fileType !=
"HEX")
5685 std::cerr<<
"readHex does not recognize filetype "
5686 <<fileType<<std::endl;
5711 int n0,n1,n2,n3,n4,n5,n6,n7,emt;
5715 meshFile>>n0>>n1>>n2>>n3>>n4>>n5>>n6>>n7>>emt;
5740 using namespace std;
5743 std::string bcFilename= std::string(filebase)+
".bc";
5744 std::ifstream bcFile(bcFilename.c_str());
5747 std::cerr<<
"readBC cannot open file "
5748 <<bcFilename<<std::endl;
5752 std::string firstWord;
5753 int eN,ebN_local,nN,flag;
5754 while (!bcFile.eof())
5757 if (firstWord ==
"FCS")
5759 bcFile>>eN>>ebN_local>>flag;
5762 if (firstWord ==
"NDS")
5776 using namespace std;
5796 using namespace std;
5815 using namespace std;
5818 int nLevelsNew = multilevelMesh.
nLevels+nLevels2add;
5819 Mesh * meshArrayTmp =
new Mesh[nLevelsNew];
5820 int** elementChildrenArrayTmp =
new int*[nLevelsNew];
5821 int** elementChildrenOffsetsTmp =
new int*[nLevelsNew];
5822 int** elementParentsArrayTmp =
new int*[nLevelsNew];
5825 for (
int i=0; i < multilevelMesh.
nLevels; i++)
5827 meshArrayTmp[i] = multilevelMesh.
meshArray[i];
5837 multilevelMesh.
meshArray = meshArrayTmp;
5843 for (
int i=multilevelMesh.
nLevels; i < nLevelsNew; i++)
5850 multilevelMesh.
nLevels = nLevelsNew;
5855 int * elementTagArray)
5857 using namespace std;
5863 int nLevelsPrev = multilevelMesh.
nLevels;
5865 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
5870 int nElements_tagged = 0;
5873 if (elementTagArray[eN] > 0)
5927 set<Node> newNodeSet;
5928 set<Node>::iterator nodeItr;
5933 if (elementTagArray[eN_parent] == 0)
5947 else if (elementTagArray[eN_parent] > 0)
5964 midpoints[0].
nN = nN_new;
5965 newNodeSet.insert(midpoints[0]);
5982 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
6005 int * elementTagArray)
6007 using namespace std;
6013 int nLevelsPrev = multilevelMesh.
nLevels;
6017 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
6018 int nElements_tagged = 0;
6021 if (elementTagArray[eN] > 0)
6027 vector<int> elementNodesArray_tmp,elementNeighborsArray_tmp,bases_tmp,
6028 elementParentsArray_tmp;
6029 vector<double> nodeArray_tmp;
6036 4*nElements_tagged));
6038 4*nElements_tagged));
6040 4*nElements_tagged);
6043 3*nElements_tagged));
6047 4*nElements_tagged);
6051 elementNodesArray_tmp.insert(elementNodesArray_tmp.begin(),
6056 elementNeighborsArray_tmp.insert(elementNeighborsArray_tmp.begin(),
6061 bases_tmp.insert(bases_tmp.begin(),
6067 nodeArray_tmp.insert(nodeArray_tmp.begin(),
6078 elementParentsArray_tmp.push_back(eN);
6080 for (
int eN_parent = 0;
6085 if (elementTagArray[eN_parent] == 1 && !refined[eN_parent])
6091 elementNodesArray_tmp,
6092 elementNeighborsArray_tmp,
6093 elementChildrenList,
6094 elementParentsArray_tmp,
6101 assert(elementParentsArray_tmp.size() ==
unsigned(nElements_new));
6104 copy(elementParentsArray_tmp.begin(),elementParentsArray_tmp.end(),
6115 if (elementChildrenList[eN].size() > 0)
6118 offset + elementChildrenList[eN].size();
6120 for (std::list<int>::iterator it = elementChildrenList[eN].begin();
6121 it != elementChildrenList[eN].end(); it++)
6156 copy(elementNodesArray_tmp.begin(),elementNodesArray_tmp.end(),
6158 copy(nodeArray_tmp.begin(),nodeArray_tmp.end(),
6160 copy(bases_tmp.begin(),bases_tmp.end(),
6185 using namespace std;
6186 int nLevels = multilevelMesh.
nLevels;
6189 cout<<
"WARNING setNewestNodeBasesToLongestEdge newestNodeBases !=NULL exiting"<<endl;
6195 cout<<
"WARNING setNewestNodeBasesToLongestEdge 2d only exiting"<<endl;
6212 int * elementTagArray)
6214 using namespace std;
6217 int nLevelsPrev = multilevelMesh.
nLevels;
6219 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
6222 set<Node> newNodeSet;
6223 set<Node>::iterator nodeItr;
6225 map<int, vector<int> > elementChildren;
6226 map<int, int> elementParents;
6228 map<int,int> elementsForBisection;
6229 map<int,vector<int> > elementsForTrisection,newElementsForTrisection,nextElementsForTrisection;
6230 set<int> elementsForUniform,newElementsForUniform,nextElementsForUniform;
6232 vector<valarray<int> > newElements;
6234 set<int> oldElements;
6235 int i = nLevelsPrev;
6239 oldElements.insert(eN_parent);
6240 if (elementTagArray[eN_parent] > 0)
6242 newElementsForUniform.insert(eN_parent);
6243 elementsForUniform.insert(eN_parent);
6248 while (!newElementsForUniform.empty() || !newElementsForTrisection.empty())
6251 for(set<int>::iterator eN_uniform_itr = newElementsForUniform.begin(); eN_uniform_itr != newElementsForUniform.end(); eN_uniform_itr++)
6253 int eN_parent = *eN_uniform_itr;
6254 oldElements.erase(eN_parent);
6255 elementChildren[eN_parent].push_back(eN_new+0);
6256 elementChildren[eN_parent].push_back(eN_new+1);
6257 elementChildren[eN_parent].push_back(eN_new+2);
6258 elementChildren[eN_parent].push_back(eN_parent);
6259 elementParents[eN_new + 0] = eN_parent;
6260 elementParents[eN_new + 1] = eN_parent;
6261 elementParents[eN_new + 2] = eN_parent;
6262 elementParents[eN_parent] = eN_parent;
6264 for(
int nN_element_0=0,nN_midpoint=0;nN_element_0<multilevelMesh.
meshArray[i-1].
nNodes_element;nN_element_0++)
6265 for(
int nN_element_1=nN_element_0+1;nN_element_1<multilevelMesh.
meshArray[i-1].
nNodes_element;nN_element_1++,nN_midpoint++)
6269 midpoints[nN_midpoint]);
6270 nodeItr = newNodeSet.find(midpoints[nN_midpoint]);
6271 if(nodeItr == newNodeSet.end())
6273 midpoints[nN_midpoint].
nN = nN_new;
6274 newNodeSet.insert(midpoints[nN_midpoint]);
6278 midpoints[nN_midpoint].
nN = nodeItr->nN;
6281 valarray<int> newElement(4);
6282 newElement[0] = eN_new;
6284 newElement[2] = midpoints[0].
nN;
6285 newElement[3] = midpoints[1].
nN;
6286 newElements.push_back(newElement);
6288 newElement[0] = eN_new;
6290 newElement[2] = midpoints[0].
nN;
6291 newElement[3] = midpoints[2].
nN;
6292 newElements.push_back(newElement);
6294 newElement[0] = eN_new;
6296 newElement[2] = midpoints[1].
nN;
6297 newElement[3] = midpoints[2].
nN;
6298 newElements.push_back(newElement);
6300 newElement[0] = eN_parent;
6301 newElement[1] = midpoints[0].
nN;
6302 newElement[2] = midpoints[1].
nN;
6303 newElement[3] = midpoints[2].
nN;
6304 newElements.push_back(newElement);
6312 ebN_neighbor_element=0;
6317 if (eN_neighbor != -1 && elementsForUniform.find(eN_neighbor) == elementsForUniform.end())
6319 if (elementsForTrisection.find(eN_neighbor) == elementsForTrisection.end())
6321 if(elementsForBisection.find(eN_neighbor) == elementsForBisection.end())
6324 int leftNode,rightNode,
6325 leftNode1,rightNode1,
6326 ebN_neighbor_element1,ebN_global1,
6327 longestEdge,longestEdge_element;
6332 longestEdge=ebN_global;
6333 longestEdge_element=ebN_neighbor_element;
6338 ebN_neighbor_element1];
6344 longestEdge = ebN_global1;
6345 longestEdge_element=ebN_neighbor_element1;
6348 if (longestEdge == ebN_global)
6350 elementsForBisection[eN_neighbor] = longestEdge_element;
6354 newElementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6355 newElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6356 elementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6357 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6362 if (elementsForBisection[eN_neighbor] != ebN_neighbor_element)
6364 newElementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6365 newElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6366 elementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6367 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6368 elementsForBisection.erase(eN_neighbor);
6374 if (elementsForTrisection[eN_neighbor][0] != ebN and elementsForTrisection[eN_neighbor][0] != ebN)
6376 nextElementsForUniform.insert(eN_neighbor);
6377 elementsForUniform.insert(eN_neighbor);
6378 elementsForTrisection.erase(eN_neighbor);
6379 newElementsForTrisection.erase(eN_neighbor);
6385 newElementsForUniform.clear();
6386 newElementsForUniform = nextElementsForUniform;
6387 nextElementsForUniform.clear();
6389 for(map<
int,vector<int> >::iterator eN_trisection_itr = newElementsForTrisection.begin(); eN_trisection_itr != newElementsForTrisection.end(); eN_trisection_itr++)
6391 int eN_parent = eN_trisection_itr->first,
6392 longestEdge_element = eN_trisection_itr->second[0],
6393 otherEdge_element = eN_trisection_itr->second[1];
6401 nodeItr = newNodeSet.find(midpoint_new);
6402 if(nodeItr == newNodeSet.end())
6404 midpoint_new.
nN = nN_new;
6405 newNodeSet.insert(midpoint_new);
6411 longestEdge_element],
6412 ebN_neighbor_element=0;
6417 if (eN_neighbor != -1 && elementsForUniform.find(eN_neighbor) == elementsForUniform.end())
6419 if (elementsForTrisection.find(eN_neighbor) == elementsForTrisection.end())
6421 if(elementsForBisection.find(eN_neighbor) == elementsForBisection.end())
6424 int leftNode,rightNode,
6425 leftNode1,rightNode1,
6426 ebN_neighbor_element1,ebN_global1,
6427 longestEdge,longestEdge_element;
6432 longestEdge=ebN_global;
6433 longestEdge_element=ebN_neighbor_element;
6438 ebN_neighbor_element1];
6445 longestEdge = ebN_global1;
6446 longestEdge_element=ebN_neighbor_element1;
6449 if (longestEdge == ebN_global)
6451 elementsForBisection[eN_neighbor] = longestEdge_element;
6455 nextElementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6456 nextElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6457 elementsForTrisection[eN_neighbor].push_back(longestEdge_element);
6458 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6463 assert(elementsForBisection[eN_neighbor] != ebN_neighbor_element);
6464 nextElementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6465 nextElementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6466 elementsForTrisection[eN_neighbor].push_back(elementsForBisection[eN_neighbor]);
6467 elementsForTrisection[eN_neighbor].push_back(ebN_neighbor_element);
6468 elementsForBisection.erase(eN_neighbor);
6473 if (elementsForTrisection[eN_neighbor][0] != longestEdge and elementsForTrisection[eN_neighbor][0] != longestEdge)
6475 newElementsForUniform.insert(eN_neighbor);
6476 elementsForUniform.insert(eN_neighbor);
6477 elementsForTrisection.erase(eN_neighbor);
6478 nextElementsForTrisection.erase(eN_neighbor);
6483 Node midpoint_other;
6487 nodeItr = newNodeSet.find(midpoint_other);
6488 assert(nodeItr != newNodeSet.end());
6490 newElementsForTrisection.clear();
6491 newElementsForTrisection = nextElementsForTrisection;
6492 nextElementsForTrisection.clear();
6495 for(map<
int,vector<int> >::iterator eN_trisection_itr = elementsForTrisection.begin(); eN_trisection_itr != elementsForTrisection.end(); eN_trisection_itr++)
6497 int eN_parent = eN_trisection_itr->first,
6498 longestEdge_element = eN_trisection_itr->second[0],
6499 otherEdge_element = eN_trisection_itr->second[1];
6507 if (otherEdge_rightNode == longestEdge_leftNode)
6509 n0 = otherEdge_rightNode;
6510 n1 = otherEdge_leftNode;
6511 n2 = longestEdge_rightNode;
6513 else if (otherEdge_rightNode == longestEdge_rightNode)
6515 n0 = otherEdge_rightNode;
6516 n1 = otherEdge_leftNode;
6517 n2 = longestEdge_leftNode;
6519 else if (otherEdge_leftNode == longestEdge_leftNode)
6521 n0 = otherEdge_leftNode;
6522 n1 = otherEdge_rightNode;
6523 n2 = longestEdge_rightNode;
6527 n0 = otherEdge_leftNode;
6528 n1 = otherEdge_rightNode;
6529 n2 = longestEdge_leftNode;
6531 oldElements.erase(eN_parent);
6532 elementChildren[eN_parent].push_back(eN_new+0);
6533 elementChildren[eN_parent].push_back(eN_new+1);
6534 elementChildren[eN_parent].push_back(eN_parent);
6535 elementParents[eN_new + 0] = eN_parent;
6536 elementParents[eN_new + 1] = eN_parent;
6537 elementParents[eN_parent] = eN_parent;
6543 nodeItr = newNodeSet.find(midpoint_new);
6544 assert(nodeItr != newNodeSet.end());
6545 midpoint_new.
nN = nodeItr->nN;
6546 Node midpoint_other;
6550 nodeItr = newNodeSet.find(midpoint_other);
6551 assert(nodeItr != newNodeSet.end());
6552 midpoint_other.
nN = nodeItr->nN;
6554 valarray<int> newElement(4);
6555 newElement[0] = eN_new;
6557 newElement[2] = midpoint_new.
nN;
6558 newElement[3] = midpoint_other.
nN;
6559 newElements.push_back(newElement);
6561 newElement[0] = eN_new;
6563 newElement[2] = midpoint_other.
nN;
6564 newElement[3] = midpoint_new.
nN;
6565 newElements.push_back(newElement);
6567 newElement[0] = eN_parent;
6569 newElement[2] = midpoint_new.
nN;
6571 newElements.push_back(newElement);
6576 for(map<int,int>::iterator eN_bisect_itr = elementsForBisection.begin(); eN_bisect_itr != elementsForBisection.end(); eN_bisect_itr++)
6578 int eN_parent = eN_bisect_itr->first;
6610 oldElements.erase(eN_parent);
6612 elementChildren[eN_parent].push_back(eN_new);
6613 elementChildren[eN_parent].push_back(eN_parent);
6614 elementParents[eN_new] = eN_parent;
6615 elementParents[eN_parent] = eN_parent;
6618 int nN_element_1,nN_element_2;
6619 bool foundMidpoint=
false;
6622 foundMidpoint=
false;
6623 nN_element_1 = (nN_element_0+1)%3;
6624 nN_element_2 = (nN_element_0+2)%3;
6628 nodeItr = newNodeSet.find(m);
6629 if(nodeItr != newNodeSet.end())
6632 valarray<int> newElement(4);
6633 newElement[0] = eN_new;
6636 newElement[3] = nodeItr->nN;
6637 newElements.push_back(newElement);
6639 newElement[0] = eN_parent;
6641 newElement[2] = nodeItr->nN;
6643 newElements.push_back(newElement);
6661 for(set<int>::iterator eN_itr=oldElements.begin();eN_itr != oldElements.end();eN_itr++)
6669 for(
vector<valarray<int> >::iterator element_itr=newElements.begin();element_itr != newElements.end();element_itr++)
6671 int eN = (*element_itr)[0];
6686 for(nodeItr=newNodeSet.begin();nodeItr!=newNodeSet.end();nodeItr++)
6703 for(
unsigned int childE=0;childE<elementChildren[eN_parent].size();childE++)
6723 int * elementTagArray)
6725 using namespace std;
6731 int nLevelsPrev = multilevelMesh.
nLevels;
6733 assert(multilevelMesh.
nLevels == nLevelsPrev+1);
6734 int nElements_tagged = 0;
6737 if (elementTagArray[eN] > 0)
6749 for (
int eN_parent = 0;
6753 if (elementTagArray[eN_parent] == 1)
6779 for (
int eN_parent = 0;
6782 int nBisectedEdges = 0;
6783 if (refined[eN_parent])
6785 for (
int ebN_local = 0;
6792 if (edgeMidNodesArray[ebN] >= 0)
6796 if (nBisectedEdges == 3)
6798 else if (nBisectedEdges == 2)
6800 else if (nBisectedEdges == 1)
6845 if (edgeMidNodesArray[ebN_parent] >= 0)
6851 midpoints[0].
nN = edgeMidNodesArray[ebN_parent];
6863 bool subdivideFailed =
false;
6864 for (
int eN_parent = 0;
6878 if (subdivideFailed)
6901 int& nElements_global,
6903 std::vector<double>& nodeArray,
6904 std::vector<int>& elementNodesArray,
6905 std::vector<int>& elementNeighborsArray,
6907 std::vector<int>& elementParentsArray,
6908 std::vector<int>& bases,
6909 std::vector<bool>& refined)
6912 bool failed =
false;
6913 const int nElementBoundaries_element = 3;
6914 const int nSpace = 2;
6921 double x[3] = {0.0,0.0,0.0};
6929 int ebN_base = bases[eN];
6930 ib[0] = ebN_base; ib[1] = (ebN_base+1)%nElementBoundaries_element; ib[2]=(ebN_base+2)%nElementBoundaries_element;
6931 for (
int i=0; i < nElementBoundaries_element; i++)
6932 IB[i] = elementNodesArray[nElementBoundaries_element*eN + ib[i]];
6935 int eN_neig = elementNeighborsArray[nElementBoundaries_element*eN + ib[0]];
6943 int newNodeNumber = nNodes_global;
6944 x[0] = 0.0; x[1] = 0.0; x[2] = 0.0;
6945 for (
int nN = 0; nN < nSpace; nN++)
6946 for (
int I = 0; I < 3; I++)
6947 x[I] += 0.5*nodeArray[3*IB[nN+1]+I];
6954 for (
int I = 0; I < 3; I++)
6955 nodeArray.push_back(x[I]);
6967 int newElementNumber = nElements_global;
6968 E1[0] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
6969 E1[1] = elementNodesArray[nElementBoundaries_element*eN + ib[1]];
6970 E1[2] = newNodeNumber;
6971 E2[0] = newNodeNumber;
6972 E2[1] = elementNodesArray[nElementBoundaries_element*eN + ib[2]];
6973 E2[2] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
6976 N1[1] = newElementNumber;
6977 N1[2] = elementNeighborsArray[nElementBoundaries_element*eN + ib[2]];
6980 N2[0] = elementNeighborsArray[nElementBoundaries_element*eN + ib[1]];
6985 assert(N2[0] < nElements_global);
6986 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
6988 if (elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] == eN)
6990 elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] =
7010 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7011 elementNodesArray[nElementBoundaries_element*eN + ebN]=E1[ebN];
7013 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7014 elementNodesArray.push_back(E2[ebN]);
7016 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7017 elementNeighborsArray[nElementBoundaries_element*eN + ebN]=N1[ebN];
7018 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7019 elementNeighborsArray.push_back(N2[ebN]);
7026 if (
unsigned(eN) < refined.size() &&
7029 elementParentsArray[eN] = eN;
7030 elementParentsArray.push_back(eN);
7031 childrenList[eN].push_back(eN);
7032 childrenList[eN].push_back(newElementNumber);
7037 assert(
unsigned(eN) < elementParentsArray.size());
7039 elementParentsArray.push_back(elementParentsArray[eN]);
7042 while (
unsigned(eN_tmp) >= refined.size())
7043 eN_tmp = elementParentsArray[eN_tmp];
7044 assert(
unsigned(elementParentsArray[eN_tmp]) < childrenList.size());
7046 childrenList[eN_tmp].push_back(newElementNumber);
7051 nElements_global += 1;
7054 if (
unsigned(eN) < refined.size())
7057 else if (elementNeighborsArray[nElementBoundaries_element*eN_neig + bases[eN_neig]] == eN)
7063 assert(eN_neig < nElements_global);
7066 int newNodeNumber = nNodes_global;
7067 x[0] = 0.0; x[1] = 0.0; x[2] = 0.0;
7068 for (
int nN = 0; nN < nSpace; nN++)
7069 for (
int I = 0; I < 3; I++)
7071 x[I] += 0.5*nodeArray[3*IB[nN+1]+I];
7074 for (
int I = 0; I < 3; I++)
7075 nodeArray.push_back(x[I]);
7094 int ebN_base_neig = bases[eN_neig];
7095 ibn[0]=ebN_base_neig; ibn[1]=(ebN_base_neig+1)%nElementBoundaries_element;
7096 ibn[2]=(ebN_base_neig+2)%nElementBoundaries_element;
7098 int newElementNumber = nElements_global;
7099 int newElementNumberNeig = nElements_global+1;
7100 E1[0] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
7101 E1[1] = elementNodesArray[nElementBoundaries_element*eN + ib[1]];
7102 E1[2] = newNodeNumber;
7103 E2[0] = newNodeNumber;
7104 E2[1] = elementNodesArray[nElementBoundaries_element*eN + ib[2]];
7105 E2[2] = elementNodesArray[nElementBoundaries_element*eN + ib[0]];
7107 N1[0] = newElementNumberNeig;
7108 N1[1] = newElementNumber;
7109 N1[2] = elementNeighborsArray[nElementBoundaries_element*eN + ib[2]];
7113 N2[0] = elementNeighborsArray[nElementBoundaries_element*eN + ib[1]];
7118 assert(N2[0] < nElements_global);
7119 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7121 if (elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] == eN)
7123 elementNeighborsArray[nElementBoundaries_element*N2[0] + ebN] =
7132 E1n[0] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[0]];
7133 E1n[1] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[1]];
7134 E1n[2] = newNodeNumber;
7135 E2n[0] = newNodeNumber;
7136 E2n[1] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[2]];
7137 E2n[2] = elementNodesArray[nElementBoundaries_element*eN_neig + ibn[0]];
7139 N1n[0] = newElementNumber;
7140 N1n[1] = newElementNumberNeig;
7141 N1n[2] = elementNeighborsArray[nElementBoundaries_element*eN_neig + ibn[2]];
7145 N2n[0] = elementNeighborsArray[nElementBoundaries_element*eN_neig + ibn[1]];
7150 assert(N2n[0] < nElements_global);
7151 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7153 if (elementNeighborsArray[nElementBoundaries_element*N2n[0] + ebN] == eN_neig)
7155 elementNeighborsArray[nElementBoundaries_element*N2n[0] + ebN] =
7156 newElementNumberNeig;
7164 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7166 elementNodesArray[nElementBoundaries_element*eN + ebN] =E1[ebN];
7167 elementNodesArray[nElementBoundaries_element*eN_neig + ebN]=E1n[ebN];
7170 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7171 elementNodesArray.push_back(E2[ebN]);
7173 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7174 elementNodesArray.push_back(E2n[ebN]);
7177 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7179 elementNeighborsArray[nElementBoundaries_element*eN + ebN] =N1[ebN];
7180 elementNeighborsArray[nElementBoundaries_element*eN_neig + ebN]=N1n[ebN];
7182 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7183 elementNeighborsArray.push_back(N2[ebN]);
7184 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7185 elementNeighborsArray.push_back(N2n[ebN]);
7194 if (
unsigned(eN) < refined.size() &&
7197 elementParentsArray[eN] = eN;
7198 elementParentsArray.push_back(eN);
7199 childrenList[eN].push_back(eN);
7200 childrenList[eN].push_back(newElementNumber);
7205 assert(
unsigned(eN) < elementParentsArray.size());
7207 elementParentsArray.push_back(elementParentsArray[eN]);
7210 while (
unsigned(eN_tmp) >= refined.size())
7211 eN_tmp = elementParentsArray[eN_tmp];
7212 assert(
unsigned(elementParentsArray[eN_tmp]) < childrenList.size());
7214 childrenList[eN_tmp].push_back(newElementNumber);
7217 if (
unsigned(eN_neig) < refined.size() &&
7220 elementParentsArray[eN_neig] = eN_neig;
7221 elementParentsArray.push_back(eN_neig);
7222 childrenList[eN_neig].push_back(eN_neig);
7223 childrenList[eN_neig].push_back(newElementNumberNeig);
7228 assert(
unsigned(eN_neig) < elementParentsArray.size());
7230 elementParentsArray.push_back(elementParentsArray[eN_neig]);
7232 int eN_tmp = eN_neig;
7233 while (
unsigned(eN_tmp) >= refined.size())
7234 eN_tmp = elementParentsArray[eN_tmp];
7235 assert(
unsigned(elementParentsArray[eN_tmp]) < childrenList.size());
7237 childrenList[eN_tmp].push_back(newElementNumberNeig);
7241 nElements_global += 2;
7244 if (
unsigned(eN) < refined.size())
7247 if (
unsigned(eN_neig) < refined.size())
7248 refined[eN_neig] =
true;
7259 elementNeighborsArray,
7261 elementParentsArray,
7269 elementNeighborsArray,
7271 elementParentsArray,
7280 std::vector<bool>& refined,
7281 std::vector<int>& edgeMidNodesArray,
7282 const int* elementNodesArray,
7283 const int* elementBoundariesArray,
7284 const int* elementNeighborsArray,
7285 const double * nodeArray)
7287 bool failed =
false;
7288 const int nElementBoundaries_element = 3;
7289 bool refinedAlready = refined[eN];
7298 eN_longest = elementBoundariesArray[eN*nElementBoundaries_element+eN_longest_local];
7299 for (
int ebN_local = 0; ebN_local < nElementBoundaries_element; ebN_local++)
7301 const int ebN = elementBoundariesArray[eN*nElementBoundaries_element+ebN_local];
7303 if (edgeMidNodesArray[ebN] >= 0 && !refinedAlready)
7310 assert(edgeMidNodesArray[ebN] < 0 || refinedAlready);
7311 if (edgeMidNodesArray[ebN] < 0)
7313 edgeMidNodesArray[ebN] = nNodes_global;
7316 int eN_neig = elementNeighborsArray[eN*nElementBoundaries_element+ebN_local];
7326 elementBoundariesArray,
7327 elementNeighborsArray,
7337 int ebN_neig,
int eN_neig,
7339 std::vector<bool>& refined,
7340 std::vector<int>& edgeMidNodesArray,
7341 const int* elementNodesArray,
7342 const int* elementBoundariesArray,
7343 const int* elementNeighborsArray,
7344 const double * nodeArray)
7350 bool failed =
false;
7352 const int nElementBoundaries_element = 3;
7369 if (edgeMidNodesArray[ebN_neig] < 0)
7371 edgeMidNodesArray[ebN_neig] = nNodes_global;
7376 assert(eN_neig >=0);
7380 int ebN_neig_longest = elementBoundariesArray[eN_neig*nElementBoundaries_element+eN_neig_longest_local];
7382 refined[eN_neig] =
true;
7383 if (edgeMidNodesArray[ebN_neig_longest] < 0)
7385 edgeMidNodesArray[ebN_neig_longest] = nNodes_global;
7388 if (ebN_longest == ebN_neig_longest)
7391 assert(ebN_longest == ebN_neig);
7395 int eN_neig_neig = elementNeighborsArray[eN_neig*nElementBoundaries_element+eN_neig_longest_local];
7397 ebN_neig_longest,eN_neig_neig,
7402 elementBoundariesArray,
7403 elementNeighborsArray,
7408 const int* elementNodesArray,
7409 const double * nodeArray)
7411 const int nElementBoundaries_element = 3;
7412 int longest = 0;
double h_longest=0.0;
7414 int nN0 = elementNodesArray[eN*nElementBoundaries_element+((ebN+1)%nElementBoundaries_element)];
7415 int nN1 = elementNodesArray[eN*nElementBoundaries_element+((ebN+2)%nElementBoundaries_element)];
7416 double len = fabs((nodeArray[3*nN1+0]-nodeArray[3*nN0+0])*(nodeArray[3*nN1+0]-nodeArray[3*nN0+0])+
7417 (nodeArray[3*nN1+1]-nodeArray[3*nN0+1])*(nodeArray[3*nN1+1]-nodeArray[3*nN0+1])+
7418 (nodeArray[3*nN1+2]-nodeArray[3*nN0+2])*(nodeArray[3*nN1+2]-nodeArray[3*nN0+2]));
7420 longest = ebN; h_longest = len;
7421 for (ebN = 1; ebN < nElementBoundaries_element; ebN++)
7423 nN0 = elementNodesArray[eN*nElementBoundaries_element+((ebN+1)%nElementBoundaries_element)];
7424 nN1 = elementNodesArray[eN*nElementBoundaries_element+((ebN+2)%nElementBoundaries_element)];
7425 len = fabs((nodeArray[3*nN1+0]-nodeArray[3*nN0+0])*(nodeArray[3*nN1+0]-nodeArray[3*nN0+0])+
7426 (nodeArray[3*nN1+1]-nodeArray[3*nN0+1])*(nodeArray[3*nN1+1]-nodeArray[3*nN0+1])+
7427 (nodeArray[3*nN1+2]-nodeArray[3*nN0+2])*(nodeArray[3*nN1+2]-nodeArray[3*nN0+2]));
7429 if (len > h_longest)
7430 { len = h_longest; longest = ebN;}
7437 int* elementParentsArray,
7438 int* elementChildrenOffsets,
7439 int* elementChildrenArray,
7440 int* elementNodesArray_child,
7441 const std::vector<int>& edgeMidNodesArray,
7442 const std::vector<bool>& refined,
7443 const int* elementNodesArray_parent,
7444 const int* elementBoundariesArray_parent,
7445 const double* nodeArray_parent)
7447 bool failed =
false,u4T=
false;
7449 const int simplexDim = 3;
7450 const int childOffset = elementChildrenOffsets[eN_parent];
7451 if (!refined[eN_parent])
7454 for (
int nN = 0; nN < simplexDim; nN++)
7456 elementNodesArray_child[eN_parent*simplexDim+nN] =
7457 elementNodesArray_parent[eN_parent*simplexDim+nN];
7460 elementParentsArray[eN_parent] = eN_parent;
7461 elementChildrenOffsets[eN_parent+1]= childOffset+1;
7462 elementChildrenArray[childOffset] = eN_parent;
7467 int nBisectedEdges = 0;
7468 int midnodes[3],vnodes[3];
7469 for (
int ebN_local = 0; ebN_local < simplexDim; ebN_local++)
7471 const int ebN = elementBoundariesArray_parent[eN_parent*simplexDim+ebN_local];
7472 vnodes[ebN_local] = elementNodesArray_parent[eN_parent*simplexDim+ebN_local];
7473 midnodes[ebN_local] = edgeMidNodesArray[ebN];
7474 if (edgeMidNodesArray[ebN] >= 0)
7477 assert(nBisectedEdges > 0);
7479 if (nBisectedEdges == 3)
7486 elementChildrenOffsets[eN_parent+1]= childOffset+4;
7490 elementNodesArray_child[eN_parent*simplexDim+0]=
7492 elementNodesArray_child[eN_parent*simplexDim+1]=
7494 elementNodesArray_child[eN_parent*simplexDim+2]=
7497 elementParentsArray[eN_parent] = eN_parent;
7498 elementChildrenArray[childOffset+0] = eN_parent;
7499 for (
int c=0;
c<simplexDim;
c++)
7501 elementNodesArray_child[eN_new*simplexDim+0]=
7503 elementNodesArray_child[eN_new*simplexDim+1]=
7504 midnodes[(
c+1)%simplexDim];
7505 elementNodesArray_child[eN_new*simplexDim+2]=
7506 midnodes[(
c+2)%simplexDim];
7508 elementParentsArray[eN_new] = eN_parent;
7509 elementChildrenArray[childOffset+
c+1] = eN_new;
7529 elementNodesArray_parent,
7531 assert(midnodes[base] >= 0);
7533 elementChildrenOffsets[eN_parent+1]= childOffset+4;
7536 elementNodesArray_child[eN_parent*simplexDim+0]=
7538 elementNodesArray_child[eN_parent*simplexDim+1]=
7539 midnodes[(base+2)%simplexDim];
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]=
7547 midnodes[(base+2)%simplexDim];
7548 elementNodesArray_child[eN_new*simplexDim+1]=
7549 vnodes[(base+1)%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]=
7559 elementNodesArray_child[eN_new*simplexDim+1]=
7560 vnodes[(base+2)%simplexDim];
7561 elementNodesArray_child[eN_new*simplexDim+2]=
7562 midnodes[(base+1)%simplexDim];
7564 elementParentsArray[eN_new] = eN_parent;
7565 elementChildrenArray[childOffset+2] = eN_new;
7568 elementNodesArray_child[eN_new*simplexDim+0]=
7569 midnodes[(base+1)%simplexDim];
7570 elementNodesArray_child[eN_new*simplexDim+1]=
7572 elementNodesArray_child[eN_new*simplexDim+2]=
7575 elementParentsArray[eN_new] = eN_parent;
7576 elementChildrenArray[childOffset+3] = eN_new;
7580 else if (nBisectedEdges == 2)
7601 elementNodesArray_parent,
7603 assert(midnodes[base] >= 0);
7605 elementChildrenOffsets[eN_parent+1]= childOffset+3;
7607 if (midnodes[(base+1)%simplexDim] < 0)
7609 assert(midnodes[(base+2)%simplexDim] >= 0);
7612 elementNodesArray_child[eN_parent*simplexDim+0]=
7613 vnodes[(base+2)%simplexDim];
7614 elementNodesArray_child[eN_parent*simplexDim+1]=
7616 elementNodesArray_child[eN_parent*simplexDim+2]=
7619 elementParentsArray[eN_parent] = eN_parent;
7620 elementChildrenArray[childOffset+0] = eN_parent;
7622 elementNodesArray_child[eN_new*simplexDim+0]=
7624 elementNodesArray_child[eN_new*simplexDim+1]=
7625 midnodes[(base+2)%simplexDim];
7626 elementNodesArray_child[eN_new*simplexDim+2]=
7629 elementParentsArray[eN_new] = eN_parent;
7630 elementChildrenArray[childOffset+1] = eN_new;
7633 elementNodesArray_child[eN_new*simplexDim+0]=
7634 midnodes[(base+2)%simplexDim];
7635 elementNodesArray_child[eN_new*simplexDim+1]=
7636 vnodes[(base+1)%simplexDim];
7637 elementNodesArray_child[eN_new*simplexDim+2]=
7640 elementParentsArray[eN_new] = eN_parent;
7641 elementChildrenArray[childOffset+2] = eN_new;
7646 assert(midnodes[(base+2)%simplexDim] < 0);
7649 elementNodesArray_child[eN_parent*simplexDim+0]=
7651 elementNodesArray_child[eN_parent*simplexDim+1]=
7652 vnodes[(base+1)%simplexDim];
7653 elementNodesArray_child[eN_parent*simplexDim+2]=
7656 elementParentsArray[eN_parent] = eN_parent;
7657 elementChildrenArray[childOffset+0] = eN_parent;
7659 elementNodesArray_child[eN_new*simplexDim+0]=
7661 elementNodesArray_child[eN_new*simplexDim+1]=
7663 elementNodesArray_child[eN_new*simplexDim+2]=
7664 midnodes[(base+1)%simplexDim];
7666 elementParentsArray[eN_new] = eN_parent;
7667 elementChildrenArray[childOffset+1] = eN_new;
7670 elementNodesArray_child[eN_new*simplexDim+0]=
7672 elementNodesArray_child[eN_new*simplexDim+1]=
7673 vnodes[(base+2)%simplexDim];
7674 elementNodesArray_child[eN_new*simplexDim+2]=
7675 midnodes[(base+1)%simplexDim];
7677 elementParentsArray[eN_new] = eN_parent;
7678 elementChildrenArray[childOffset+2] = eN_new;
7682 else if (nBisectedEdges == 1)
7694 elementNodesArray_parent,
7697 elementChildrenOffsets[eN_parent+1]= childOffset+2;
7699 assert(midnodes[base] >= 0);
7701 elementNodesArray_child[eN_parent*simplexDim+0]=
7703 elementNodesArray_child[eN_parent*simplexDim+1]=
7704 vnodes[(base+1)%simplexDim];
7705 elementNodesArray_child[eN_parent*simplexDim+2]=
7708 elementParentsArray[eN_parent] = eN_parent;
7709 elementChildrenArray[childOffset+0] = eN_parent;
7711 elementNodesArray_child[eN_new*simplexDim+0]=
7713 elementNodesArray_child[eN_new*simplexDim+1]=
7714 vnodes[(base+2)%simplexDim];
7715 elementNodesArray_child[eN_new*simplexDim+2]=
7718 elementParentsArray[eN_new] = eN_parent;
7719 elementChildrenArray[childOffset+1] = eN_new;
7724 int child_start = elementChildrenOffsets[eN_parent];
7725 int child_end = elementChildrenOffsets[eN_parent+1];
7742 #ifdef MWF_HACK_2D_COARSE
7743 int findGlueNeighbor2d(
int eN,
7744 std::vector<int>& elementNeighborsArray,
7745 std::vector<int>& elementParentsArray)
7750 const int nElementBoundaries_element = 3;
7752 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7754 const int eN_neig = elementNeighborsArray[eN*nElementBoundaries_element+ebN];
7755 if (eN_neig >= 0 && elementParentsArray[eN_neig] == elementParentsArray[eN])
7763 int findT2Neighbor(
int eN,
int eN_base,
7764 std::vector<int>& elementNeighborsArray,
7765 std::vector<int>& elementParentsArray)
7774 for (
int ebN = 0; ebN < nElementBoundaries_element; ebN++)
7776 const int eN_neig = elementNeighborsArray[eN*nElementBoundaries_element+ebN];
7777 if (ebN != eN_base && eN_neig >= 0 &&
7778 elementParentsArray[eN_neig] != elementParentsArray[eN])
7786 bool newestNodeGlue(
int eN,
7787 std::vector<bool>& mayCoarsen,
7788 int& nElements_global,
7790 std::vector<double>& nodeArray,
7791 std::vector<int>& elementNodesArray,
7792 std::vector<int>& elementNeighborsArray,
7793 std::vector<int>& elementParentsArray,
7794 std::vector<int>& bases,
7795 std::vector<bool>& coarsened)
7797 using namespace std;
7798 bool failed =
false;
7799 if (!mayCoarsen[eN])
7800 {failed =
true;
return failed;}
7802 const int nElementBoundaries_element = 3;
7803 bool connectable =
false;
7804 int eN_base = bases[eN];
7805 int nN_global_base = elementNodesArray[eN*nElementBoundaries_element+eN_base];
7806 int eN_glue = -1, eN_t2 = -1;
7807 eN_glue = findGlueNeighbor(eN,
7808 elementNeighborsArray,
7809 elementParentsArray);
7810 eN_t2 = findT2Neighbor(eN,eN_base,
7811 elementNeighborsArray,
7812 elementParentsArray);
7814 {failed =
true;
return failed;}
7815 if (eN_t2 >= 0 && elementNodesArray[eN_t2*nElementBoundaries_element+bases[eN_t2]] != nN_global_base)
7816 failed = newestNodeGlue(eN_t2,
7822 elementNeighborsArray,
7823 elementParentsArray,
7829 connectable = (nN_global_base == elementNodesArray[eN_glue*nElementBoundaries_element+bases[eN_glue]]);
7831 failed = newestNodeGlue(eN_glue,
7837 elementNeighborsArray,
7838 elementParentsArray,
7844 eN_glue = findGlueNeighbor(eN,
7845 elementNeighborsArray,
7846 elementParentsArray);
7848 { failed =
true;
return failed;}
7849 eN_t2 = findT2Neighbor(eN,eN_base,
7850 elementNeighborsArray,
7851 elementParentsArray);
7852 bool t2connectable =
false;
7855 int eN_t2_glue = -1;
7857 eN_t2_glue = findGlueNeighbor(eN_t2,
7858 elementNeighborsArray,
7859 elementParentsArray);
7861 if (eN_t2_glue < 0 && eN_t2 >= 0)
7862 {failed =
true;
return failed;}
7863 t2connectable = (eN_t2 >= 0 &&
7864 (elementNodesArray[eN_t2*nElementBoundaries_element+bases[eN_t2]] ==
7865 elementNodesArray[eN_t2_glue*nElementBoundaries_element+bases[eN_t2_glue]]));
7866 if (eN_t2 >= 0 && !t2connectable)
7867 failed = newestNodeGlue(eN_t2_glue,
7873 elementNeighborsArray,
7874 elementParentsArray,
7879 t2connectable =
false;
7881 eN_t2_glue = findGlueNeighbor(eN_t2,
7882 elementNeighborsArray,
7883 elementParentsArray);
7884 if (eN_t2 >= 0 && eN_t2_glue >= 0)
7885 t2connectable = (elementNodesArray[eN_t2*nElementBoundaries_element+bases[eN_t2]] ==
7886 elementNodesArray[eN_t2_glue*nElementBoundaries_element+bases[eN_t2_glue]]);
7888 failed = glueElements(eN_t2,eN_t2_glue,
7894 elementNeighborsArray,
7895 elementParentsArray,
7900 failed = glueElements(eN,eN_glue,
7906 elementNeighborsArray,
7907 elementParentsArray,
7914 bool glueElements(
int eN,
int eN_glue,
7915 std::vector<bool>& mayCoarsen,
7916 int& nElements_global,
7918 std::vector<double>& nodeArray,
7919 std::vector<int>& elementNodesArray,
7920 std::vector<int>& elementNeighborsArray,
7921 std::vector<int>& elementParentsArray,
7922 std::vector<int>& bases,
7923 std::vector<bool>& coarsened)
7927 const int nElementBoundaries_element = 3;
7928 assert(eN >= 0 && eN_glue >= 0);
7929 if (!mayCoarsen[eN] || !mayCoarsen[eN_glue])
7931 bool connectable = (elementNodesArray[eN*nElementBoundaries_element+bases[eN]] ==
7932 elementNodesArray[eN_glue*nElementBoundaries_element+bases[eN_glue]]);
7947 base_parent = (bases[eN]+1) % nElementBoundaries_element;
7948 int nN_global_gone = elementNodesArray[eN*nElementBoundaries_element+bases[eN]];
7950 for (
int nN = 0; nN < nElementBoundaries_element; nN++)
7952 if (nN == bases[eN])
7953 elementNodesArray[eN*nElementBoundaries_element+ebN] = (bases[eN_glue]+1) % nElementBoundaries_element;
7955 coarsened[eN] =
true;
7958 #endif //2d HACK COARSEN