29 std::array<MFloat, nDim> res) {
35 for(
MInt v = 0; v < n; v++) {
36 for(
MInt i = 0; i < nDim_; i++) {
37 res[i] += (*vertices)[indices[v]].coordinates[i];
42 for(
MInt i = 0; i < nDim_; i++)
53 const MInt noPoints = (signed)(*face).vertices.size();
54 array<MFloat, 3> faceMidPoint{};
56 for(
MInt i = 0; i < noPoints; i++) {
57 const MInt vertex = (*face).vertices[i];
58 for(
MInt j = 0; j < nDim; j++) {
59 faceMidPoint[j] += (*vertices)[vertex].coordinates[j];
63 for(
MInt j = 0; j < nDim; j++) {
64 faceMidPoint[j] /= noPoints;
67 array<MFloat, 3> center{};
72 center = faceMidPoint;
74 const MInt vertex0 = (*face).vertices[0];
75 for(
MInt i = 1; i < noPoints - 1; i++) {
76 const MInt vertex1 = (*face).vertices[i];
77 const MInt vertex2 = (*face).vertices[i + 1];
82 for(
MInt j = 0; j < nDim; j++) {
83 a[j] = (*vertices)[vertex1].coordinates[j] - (*vertices)[vertex0].coordinates[j];
84 b[j] = (*vertices)[vertex2].coordinates[j] - (*vertices)[vertex0].coordinates[j];
85 c[j] = MC[j] - (*vertices)[vertex0].coordinates[j];
86 centerTri[j] = (*vertices)[vertex0].coordinates[j] + (*vertices)[vertex1].coordinates[j]
87 + (*vertices)[vertex2].coordinates[j];
88 centerTetra[j] = centerTri[j] + MC[j];
91 aCrossB[0] =
a[1] *
b[2] -
a[2] *
b[1];
92 aCrossB[1] =
a[2] *
b[0] -
a[0] *
b[2];
93 aCrossB[2] =
a[0] *
b[1] -
a[1] *
b[0];
96 for(
MInt j = 0; j < nDim; j++) {
97 aCrossBTimesC += aCrossB[j] * c[j];
98 normACrossB += aCrossB[j] * aCrossB[j];
100 centerTetra[j] *= F1B4;
102 normACrossB = sqrt(normACrossB);
103 MFloat volumeTetra = F1B6 * aCrossBTimesC;
104 MFloat areaTri = F1B2 * normACrossB;
108 for(
MInt j = 0; j < nDim; j++) {
109 center[j] += areaTri * centerTri[j];
110 XC[j] += volumeTetra * centerTetra[j];
114 for(
MInt j = 0; j < nDim; j++) {
119 center = faceMidPoint;
123 ASSERT(!(std::isnan(area)),
"");
125 for(
MInt j = 0; j < nDim; j++) {
126 face->
center[j] = center[j];
127 ASSERT(!(std::isnan(center[j])),
"");
135 std::vector<polyFace>* faces,
136 const std::vector<polyVertex>* vertices) {
137 for(
MInt cC = 0; (unsigned)cC < cutCells->size(); cC++) {
139 MFloat center[3] = {F0, F0, F0};
141 MFloat M[3] = {F0, F0, F0};
143 for(
MInt i = 0; (unsigned)i < cutCell->faces.size(); i++) {
144 const MInt face = cutCell->faces[i];
145 for(
MInt j = 0; j < (signed)(*faces)[face].vertices.size(); j++) {
146 const MInt vertex = (*faces)[face].vertices[j];
147 for(
MInt k = 0; k < nDim; k++) {
148 M[k] += (*vertices)[vertex].coordinates[k];
154 for(
MInt k = 0; k < nDim; k++) {
158 for(
MInt i = 0; (unsigned)i < cutCell->faces.size(); i++) {
159 polyFace* face = &(*faces)[cutCell->faces[i]];
160 compFaceIntegrals_pyraBased3(face, vertices, M, &volume, center);
163 cutCell->volume = volume;
164 if(fabs(volume) > 1e2 * m_eps) {
165 for(
MInt i = 0; i < nDim; i++) {
166 cutCell->center[i] = center[i] / volume;
169 cutCell->volume = fabs(volume);
170 for(
MInt i = 0; i < nDim; i++) {
171 cutCell->center[i] = M[i];
182 const std::vector<polyVertex>* vertices,
MInt set) {
183 const MChar* fileName =
"cell_";
184 stringstream fileName2;
185 fileName2 << fileName << grid().tree().globalId(cellId) <<
"_s" << set <<
".vtk";
187 ofl.open((fileName2.str()).c_str(), ofstream::trunc);
191 ofl.setf(ios::fixed);
194 ofl <<
"# vtk DataFile Version 3.0" << endl
195 <<
"MAIAD cutsurface file" << endl
198 <<
"DATASET UNSTRUCTURED_GRID" << endl
201 ofl <<
"POINTS " << (*vertices).size() <<
" float" << endl;
203 for(
MInt v = 0; (unsigned)v < (*vertices).size(); v++) {
204 for(
MInt i = 0; i < 3; i++)
205 ofl << (*vertices)[v].coordinates[i] <<
" ";
212 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
213 numPoints += (*faces)[f].vertices.size();
216 ofl <<
"CELLS " << noFaces <<
" " << noFaces + numPoints << endl;
218 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
219 MInt noEdges = (signed)(*faces)[f].vertices.size();
220 ofl << noEdges <<
" ";
222 for(
MInt i = noEdges - 1; i >= 0; i--) {
223 ofl << (*faces)[f].vertices[i] <<
" ";
229 ofl <<
"CELL_TYPES " << noFaces << endl;
230 for(
MInt i = 0; i < noFaces; i++) {
234 ofl <<
"CELL_DATA " << noFaces << endl;
235 ofl <<
"SCALARS bodyId int 1" << endl;
236 ofl <<
"LOOKUP_TABLE default" << endl;
237 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
238 ofl << (*faces)[f].bodyId << endl;
241 ofl <<
"SCALARS face int 1" << endl;
242 ofl <<
"LOOKUP_TABLE default" << endl;
243 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
247 ofl <<
"SCALARS faceId int 1" << endl;
248 ofl <<
"LOOKUP_TABLE default" << endl;
249 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
250 ofl << (*faces)[f].faceId << endl;
253 ofl <<
"SCALARS faceType int 1" << endl;
254 ofl <<
"LOOKUP_TABLE default" << endl;
255 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
256 ofl << (*faces)[f].faceType << endl;
259 ofl <<
"SCALARS cutCell int 1" << endl;
260 ofl <<
"LOOKUP_TABLE default" << endl;
261 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
262 ofl << (*faces)[f].cutCell << endl;
265 ofl <<
"SCALARS isLine int 1" << endl;
266 ofl <<
"LOOKUP_TABLE default" << endl;
267 for(
MInt f = 0; (unsigned)f < (*faces).size(); f++) {
268 ofl << (*faces)[f].isLine << endl;
271 ofl <<
"POINT_DATA " << (*vertices).size() << endl;
272 ofl <<
"SCALARS point int 1" << endl;
273 ofl <<
"LOOKUP_TABLE default" << endl;
274 for(
MInt v = 0; (unsigned)v < (*vertices).size(); v++) {
278 ofl <<
"SCALARS pointId int 1" << endl;
279 ofl <<
"LOOKUP_TABLE default" << endl;
280 for(
MInt v = 0; (unsigned)v < (*vertices).size(); v++) {
281 ofl << (*vertices)[v].pointId << endl;
284 ofl <<
"SCALARS pointType int 1" << endl;
285 ofl <<
"LOOKUP_TABLE default" << endl;
286 for(
MInt v = 0; (unsigned)v < (*vertices).size(); v++) {
287 ofl << (*vertices)[v].pointType << endl;
297 MInt cellId = cutCellData[cutc].cellId;
299 const MChar* fileName =
"cell-Info_";
300 stringstream fileName2;
301 fileName2 << fileName << grid().tree().globalId(cellId) <<
"_Id" << nameId <<
".txt";
304 ofl.open((fileName2.str()).c_str(), ofstream::trunc);
306 const MInt lastDim = nDim - 1;
309 ofl.setf(ios::fixed);
312 ofl <<
"Writing Cell " << grid().tree().globalId(cellId) <<
" " << cellId <<
" " << cutc <<
" vol "
313 << cutCellData[cutc].volume <<
" SplitChilds " << cutCellData[cutc].noSplitChilds <<
" "
314 << cutCellData[cutc].splitParentId << endl
316 for(
MInt dir = 0; dir < m_noDirs; dir++) {
317 ofl <<
"Dir " << dir <<
" external " << cutCellData[cutc].externalFaces[dir] <<
" faces "
318 << cutCellData[cutc].noFacesPerCartesianDir[dir] << endl
322 ofl <<
" Cartesian-Surf: " << cutCellData[cutc].noCartesianSurfaces << endl << endl;
324 for(
MInt id = 0;
id < cutCellData[cutc].noCartesianSurfaces;
id++) {
325 ofl <<
"Id " <<
id <<
" dir " << cutCellData[cutc].cartFaceDir[
id] <<
" area "
326 << cutCellData[cutc].cartFaceArea[
id] <<
" coord " << cutCellData[cutc].cartFaceCentroid[
id][0]
327 << cutCellData[cutc].cartFaceCentroid[
id][1] << cutCellData[cutc].cartFaceCentroid[
id][lastDim] << endl
331 ofl <<
" Bndry-Surfaces: " << cutCellData[cutc].noBoundarySurfaces << endl << endl;
333 for(
MInt id = 0;
id < cutCellData[cutc].noBoundarySurfaces;
id++) {
334 ofl <<
"Id " <<
id <<
" normal " << cutCellData[cutc].boundarySurfaceNormal[
id][0]
335 << cutCellData[cutc].boundarySurfaceNormal[
id][1] << cutCellData[cutc].boundarySurfaceNormal[
id][lastDim]
336 <<
" center " << cutCellData[cutc].boundarySurfaceCentroid[
id][0]
337 << cutCellData[cutc].boundarySurfaceCentroid[
id][1] << cutCellData[cutc].boundarySurfaceCentroid[
id][lastDim]
338 <<
" area " << cutCellData[cutc].boundarySurfaceArea[
id] << endl
353 static constexpr MInt edgeCode[24] = {8, 10, 0, 4, 9, 11, 1, 5, 2, 6, 8, 9, 3, 7, 10, 11, 0, 1, 2, 3, 4, 5, 6, 7};
354 static constexpr MInt faceCode[24] = {0, 4, 1, 4, 2, 4, 3, 4, 0, 5, 1, 5, 2, 5, 3, 5, 0, 2, 1, 2, 0, 3, 1, 3};
355 static constexpr MInt pointCode[6][4] = {{0, 2, 6, 4}, {1, 3, 7, 5}, {0, 1, 5, 4},
356 {2, 3, 7, 6}, {0, 1, 3, 2}, {4, 5, 7, 6}};
357 static constexpr MInt edgeCode2[6][4] = {{0, 10, 4, 8}, {1, 11, 5, 9}, {2, 9, 6, 8},
358 {3, 11, 7, 10}, {2, 1, 3, 0}, {6, 5, 7, 4}};
359 static constexpr MFloat signStencil[8][3] = {{-F1, -F1, -F1}, {F1, -F1, -F1}, {-F1, F1, -F1}, {F1, F1, -F1},
360 {-F1, -F1, F1}, {F1, -F1, F1}, {-F1, F1, F1}, {F1, F1, F1}};
362 static constexpr MBool improvedCentroidComputation =
false;
366 MBool cutFaces[m_noDirs];
369 MBool shapeIsInside[m_noDirs];
374 const MFloat epsNormal = 1e-15;
378 MInt face0Space, face0Side;
382 MInt spaceId, spaceId1, spaceId2;
383 MInt space0, space1, space0Relative, space1Relative;
386 MInt noNonFluidEdges[m_noDirs];
391 MFloat cellHalfLength[nDim];
394 MFloat faceVolume[m_noDirs];
395 MFloat gridFaceVolume[m_noDirs];
401 MFloat minX, maxX, minY, maxY;
402 MFloat negativePoint[3] = {0, 0, 0};
403 MFloat oppositePoint[3] = {0, 0, 0};
404 MFloat point[3] = {0, 0, 0};
405 MFloat pointP[3] = {0, 0, 0};
406 MFloat pointM[3] = {0, 0, 0};
407 MFloat pointW[3] = {0, 0, 0};
408 MFloat pointE[3] = {0, 0, 0};
409 MFloat rectVolume, triVolume;
410 MFloat trianglePoint[3] = {0, 0, 0};
411 MFloat trianglePointP[3] = {0, 0, 0};
412 MFloat trianglePointM[3] = {0, 0, 0};
413 MFloat trianglePointW[3] = {0, 0, 0};
414 MFloat trianglePointE[3] = {0, 0, 0};
415 MFloat cutLineCentroid[6][3];
417 MFloat faceCentroid[6][3];
418 MFloat triangleCentroid[6][3];
419 MFloat normalVector[6][3];
424 for(
MInt f = 0; f < m_noDirs; f++) {
436 for(
MInt i = 0; i < noEdges; i++) {
437 cutEdgesAll[i] =
false;
438 nonFluidSideEdge[i] =
false;
439 cutPointOnEdge[i] = -1;
441 for(
MInt i = 0; i < m_noDirs; i++) {
443 shapeIsInside[i] =
false;
447 for(
MInt i = 0; i < nDim; i++) {
448 cellLength[i] = a_cellLengthAtCell(cellId);
449 cellHalfLength[i] = F1B2 * cellLength[i];
455 cutPointOnEdge[cutCell.
cutEdges[cp]] = cp;
466 for(
MInt face = 0; face < m_noDirs; face++) {
469 for(
MInt e = 0; e < 4; e++) {
470 MInt cp = cutPointOnEdge[edgeCode[4 * face + e]];
474 cutEdgesAll[edgeCode[4 * face + e]] =
true;
480 cutFaces[face] =
true;
483 MBool faceIsActive =
true;
484 if(n == 1 || n > 2) {
485 cerr <<
"** ERROR create cut face: " << n <<
" cut points for face " << face <<
" body "
487 stringstream errorMessage;
488 errorMessage <<
"** ERROR create cut face: " << n <<
" cut points for face " << face << endl;
489 cerr <<
"cell " <<
cellId <<
" cog " << a_coordinate(cellId, 0) <<
" " << a_coordinate(cellId, 1) <<
" "
490 << a_coordinate(cellId, 2) << endl;
495 spaceId1 = (spaceId + 1) % nDim;
496 spaceId2 = (spaceId1 + 1) % nDim;
498 gridFaceVolume[face] = cellLength[spaceId1] * cellLength[spaceId2];
499 faceVolume[face] = gridFaceVolume[face];
502 for(
MInt i = 0; i < nDim; i++) {
503 point[i] = a_coordinate(cellId, i) - cellHalfLength[i];
504 negativePoint[i] = point[i];
505 oppositePoint[i] = point[i] + cellLength[i];
509 point[spaceId] += (
MFloat)sideId * cellLength[spaceId];
510 oppositePoint[spaceId] = point[spaceId];
513 for(
MInt i = 0; i < nDim; i++) {
514 pointP[i] = point[i];
515 pointM[i] = point[i];
516 pointW[i] = point[i];
517 pointE[i] = point[i];
519 pointP[spaceId1] += eps;
520 pointM[spaceId1] -= eps;
521 pointW[spaceId2] += eps;
522 pointE[spaceId2] -= eps;
525 if(sideId) corner =
IPOW2(spaceId);
529 MInt cornerBR = corner +
IPOW2(spaceId1);
530 MInt cornerTL = corner +
IPOW2(spaceId2);
537 const MInt pointIds[4] = {corner, cornerBR, cornerOP, cornerTL};
543 faceIsActive =
false;
547 faceCentroid[face][spaceId] = point[spaceId];
551 for(
MInt i = 0; i < nDim; i++) {
552 delta[face][i] = cutCell.
cutPoints[cutPoint[0]][i] - cutCell.
cutPoints[cutPoint[1]][i];
554 positiveSlope = delta[face][spaceId1] * delta[face][spaceId2] > 0;
556 for(
MInt i = 0; i < nDim; i++) {
557 delta[face][i] = fabs(delta[face][i]);
561 cutLineLength[face] = sqrt(
POW2(delta[face][spaceId1]) +
POW2(delta[face][spaceId2]));
564 for(
MInt i = 0; i < nDim; i++) {
565 cutLineCentroid[face][i] = F1B2 * (cutCell.
cutPoints[cutPoint[0]][i] + cutCell.
cutPoints[cutPoint[1]][i]);
573 space0Relative = cutEdges[0] / 2;
574 space1Relative = cutEdges[1] / 2;
575 space0 = (space0Relative + spaceId + 1) % nDim;
576 side0 = cutEdges[0] % 2;
577 space1 = (space1Relative + spaceId + 1) % nDim;
578 side1 = cutEdges[1] % 2;
581 normalVector[face][spaceId] = 0;
582 normalVector[face][spaceId1] = delta[face][spaceId2] / cutLineLength[face];
583 normalVector[face][spaceId2] = delta[face][spaceId1] / cutLineLength[face];
590 if(space0 != space1) {
592 faceVolume[face] = F1B2 * (delta[face][spaceId1] * delta[face][spaceId2]);
594 trianglePoint[space0] = point[space0] + cellLength[space0] * (
MFloat)side0;
595 trianglePoint[space1] = point[space1] + cellLength[space1] * (
MFloat)side1;
596 trianglePoint[spaceId] = point[spaceId];
598 for(
MInt i = 0; i < nDim; i++) {
599 trianglePointP[i] = trianglePoint[i];
600 trianglePointM[i] = trianglePoint[i];
601 trianglePointW[i] = trianglePoint[i];
602 trianglePointE[i] = trianglePoint[i];
604 trianglePointP[spaceId1] += eps;
605 trianglePointM[spaceId1] -= eps;
606 trianglePointW[spaceId2] += eps;
607 trianglePointE[spaceId2] -= eps;
609 faceCentroid[face][space0] = trianglePoint[space0] + F2B3 * (F1B2 - (
MFloat)side0) * delta[face][space0];
610 faceCentroid[face][space1] = trianglePoint[space1] + F2B3 * (F1B2 - (
MFloat)side1) * delta[face][space1];
612 for(
MInt i = 0; i < nDim; i++) {
613 triangleCentroid[face][i] = faceCentroid[face][i];
617 MInt triangleCorner = 0;
618 triangleCorner += sideId *
IPOW2(spaceId);
619 triangleCorner += side0 *
IPOW2(space0);
620 triangleCorner += side1 *
IPOW2(space1);
627 shapeIsInside[face] =
true;
628 cutOutVolume = faceVolume[face];
629 faceVolume[face] = gridFaceVolume[face] - faceVolume[face];
630 faceCentroid[face][space0] =
631 F1 / faceVolume[face]
632 * ((gridFaceVolume[face] * a_coordinate(cellId, space0)) - (cutOutVolume * faceCentroid[face][space0]));
633 faceCentroid[face][space1] =
634 F1 / faceVolume[face]
635 * ((gridFaceVolume[face] * a_coordinate(cellId, space1)) - (cutOutVolume * faceCentroid[face][space1]));
638 noNonFluidEdges[face] = 0;
640 for(
MInt e = 0; e < 4; e++) {
641 if(nonFluidSideEdge[edgeCode[4 * face + e]]) {
642 cerr <<
"Error non-fluid edge was true for edge " << edgeCode[4 * face + e] <<
" face " << face << endl;
643 cerr <<
"All face data " << noNonFluidEdges[0] << noNonFluidEdges[1] << noNonFluidEdges[2]
644 << noNonFluidEdges[3] << noNonFluidEdges[4] << noNonFluidEdges[5] << endl;
645 cerr <<
"All edge data " << nonFluidSideEdge[edgeCode[0]] << nonFluidSideEdge[edgeCode[1]]
646 << nonFluidSideEdge[edgeCode[2]] << nonFluidSideEdge[edgeCode[3]] <<
" "
647 << nonFluidSideEdge[edgeCode[4]] << nonFluidSideEdge[edgeCode[5]] << nonFluidSideEdge[edgeCode[6]]
648 << nonFluidSideEdge[edgeCode[7]] <<
" " << nonFluidSideEdge[edgeCode[8]]
649 << nonFluidSideEdge[edgeCode[9]] << nonFluidSideEdge[edgeCode[10]]
650 << nonFluidSideEdge[edgeCode[11]] << endl;
651 cerr <<
"cell " <<
cellId <<
" " << a_coordinate(cellId, 0) <<
" " << a_coordinate(cellId, 1) <<
" "
652 << a_coordinate(cellId, 2) << endl;
655 nonFluidSideEdge[edgeCode[4 * face + e]] =
false;
661 noNonFluidEdges[face] = 2;
662 nonFluidSideEdge[edgeCode[4 * face + 2 * space0Relative + ((side0 + 1) % 2)]] =
true;
663 nonFluidSideEdge[edgeCode[4 * face + 2 * space1Relative + ((side1 + 1) % 2)]] =
true;
671 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
672 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
675 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
680 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
683 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
684 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
693 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
698 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
706 noNonFluidEdges[face] = 1;
707 if(space0 == spaceId1) {
709 faceVolume[face] = cellLength[spaceId1] * (cutLineCentroid[face][spaceId2] - point[spaceId2]);
713 faceVolume[face] = gridFaceVolume[face] - faceVolume[face];
714 nonFluidSideEdge[edgeCode[4 * face + 2]] =
true;
719 maxY = cutLineCentroid[face][spaceId2] + F1B2 * delta[face][spaceId2];
720 cc0[0] = a_coordinate(cellId, spaceId1);
721 cc0[1] = F1B2 * (oppositePoint[spaceId2] + maxY);
723 cc1[0] = point[spaceId1] + F1B3 * cellLength[spaceId1];
724 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
726 cc1[0] = point[spaceId1] + F2B3 * cellLength[spaceId1];
728 cc1[1] = maxY - F1B3 * delta[face][spaceId2];
729 rectVolume = (oppositePoint[spaceId2] - maxY) * cellLength[spaceId1];
730 triVolume = faceVolume[face] - rectVolume;
732 nonFluidSideEdge[edgeCode[4 * face + 3]] =
true;
737 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
738 minY = cutLineCentroid[face][spaceId2] - F1B2 * delta[face][spaceId2];
739 cc0[0] = a_coordinate(cellId, spaceId1);
740 cc0[1] = F1B2 * (point[spaceId2] + minY);
742 cc1[0] = point[spaceId1] + F2B3 * cellLength[spaceId1];
744 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
745 cc1[0] = point[spaceId1] + F1B3 * cellLength[spaceId1];
747 cc1[1] = minY + F1B3 * delta[face][spaceId2];
748 rectVolume = (minY - point[spaceId2]) * cellLength[spaceId1];
749 triVolume = faceVolume[face] - rectVolume;
753 faceVolume[face] = cellLength[spaceId2] * (cutLineCentroid[face][spaceId1] - point[spaceId1]);
757 faceVolume[face] = gridFaceVolume[face] - faceVolume[face];
758 nonFluidSideEdge[edgeCode[4 * face]] =
true;
763 maxX = cutLineCentroid[face][spaceId1] + F1B2 * delta[face][spaceId1];
764 cc0[0] = F1B2 * (oppositePoint[spaceId1] + maxX);
765 cc0[1] = a_coordinate(cellId, spaceId2);
766 cc1[0] = maxX - F1B3 * delta[face][spaceId1];
768 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
769 cc1[1] = point[spaceId2] + F1B3 * cellLength[spaceId2];
771 cc1[1] = point[spaceId2] + F2B3 * cellLength[spaceId2];
773 rectVolume = (oppositePoint[spaceId1] - maxX) * cellLength[spaceId2];
774 triVolume = faceVolume[face] - rectVolume;
776 nonFluidSideEdge[edgeCode[4 * face + 1]] =
true;
781 normalVector[face][spaceId1] = -normalVector[face][spaceId1];
782 minX = cutLineCentroid[face][spaceId1] - F1B2 * delta[face][spaceId1];
783 cc0[0] = F1B2 * (point[spaceId1] + minX);
784 cc0[1] = a_coordinate(cellId, spaceId2);
785 cc1[0] = minX + F1B3 * delta[face][spaceId1];
787 cc1[1] = point[spaceId2] + F2B3 * cellLength[spaceId2];
789 normalVector[face][spaceId2] = -normalVector[face][spaceId2];
790 cc1[1] = point[spaceId2] + F1B3 * cellLength[spaceId2];
792 rectVolume = (minX - point[spaceId1]) * cellLength[spaceId2];
793 triVolume = faceVolume[face] - rectVolume;
797 faceCentroid[face][spaceId1] = (rectVolume * cc0[0] + triVolume * cc1[0]) / faceVolume[face];
798 faceCentroid[face][spaceId2] = (rectVolume * cc0[1] + triVolume * cc1[1]) / faceVolume[face];
801 for(
MInt i = 0; i < nDim; i++) {
812 for(
MInt e = 0; e < 4; e++) {
813 nonFluidSideEdge[edgeCode[4 * face + e]] =
true;
818 for(
MInt i = 0; i < nDim; i++) {
840 for(
MInt face = 0; face < m_noDirs; face++) {
841 MBool nonFluidSide =
true;
842 for(
MInt e = 0; e < 4; e++) {
843 if(!nonFluidSideEdge[edgeCode[4 * face + e]]) {
844 nonFluidSide =
false;
851 faceVolume[face] = F0;
861 cerr <<
"case 0 " << endl;
862 cerr << a_coordinate(cellId, 0) <<
" " << a_coordinate(cellId, 1) <<
" " << a_coordinate(cellId, 2) <<
" "
863 << cellHalfLength[0] <<
" " << grid().tree().level(cellId) << endl;
875 for(
MInt face = 0; face < m_noDirs; face++) {
888 for(
MInt e = 0; e < 4; e++) {
890 if(cutEdgesAll[edgeCode[4 * face0 + e]]) {
891 cutEdges[edgeCounter] = e;
897 if(edgeCounter != 2) cerr <<
"ERROR create cut face, tetraeter shape, not enough cut edges" << endl;
899 for(
MInt e = 0; e < 4; e++) {
900 if(cutEdgesAll[edgeCode[4 * face1 + e]]) {
901 if(faceCode[2 * edgeCode[4 * face1 + e]] != face0 && faceCode[2 * edgeCode[4 * face1 + e] + 1] != face0) {
910 if(cutCell.
cutEdges[cp] == edgeCode[4 * face1 + cutEdges[2]]) {
915 face0Space = face0 / 2;
916 face0Side = face0 % 2;
920 for(
MInt i = 0; i < nDim; i++) {
925 h = F2 * (F1B2 - (
MFloat)face0Side)
926 * (cutCell.
cutPoints[cutPoint[0]][face0Space]
927 - (negativePoint[face0Space] + (
MFloat)face0Side * cellLength[face0Space]));
928 if(shapeIsInside[face0]) {
929 cutCell.
volume = F1B3 * (gridFaceVolume[face0] - faceVolume[face0]) * h;
931 cutCell.
volume = F1B3 * faceVolume[face0] * h;
935 if(shapeIsInside[face0]) {
936 cutOutVolume = cutCell.
volume;
937 cutCell.
volume = a_gridCellVolume(cellId) - cutCell.
volume;
938 for(
MInt i = 0; i < nDim; i++) {
941 * (a_gridCellVolume(cellId) * a_coordinate(cellId, i) - cutOutVolume * cutCell.
volumetricCentroid[i]);
945 for(
MInt i = 0; i < nDim; i++) {
950 for(
MInt i = 0; i < nDim; i++) {
963 for(
MInt i = 0; i < nDim; i++) {
968 for(
MInt i = 0; i < nDim; i++) {
970 F2B3 * cutLineCentroid[face0][i] + F1B3 * cutCell.
cutPoints[cutPoint[0]][i];
980 cutLineLength[6] = F0;
981 for(
MInt face = 0; face < m_noDirs; face += 2) {
983 if(cutFaces[face + 1]) {
984 cutLineLength[7] = cutLineLength[face] + cutLineLength[face + 1];
985 if(cutLineLength[7] > cutLineLength[6]) {
988 cutLineLength[6] = cutLineLength[7];
993 if(face0 < 0 || face1 < 0) {
994 cerr <<
"Error: opposite cut faces not found" << endl;
999 spaceId1 = (spaceId + 1) % nDim;
1000 spaceId2 = (spaceId1 + 1) % nDim;
1004 Fvol = F1 / (faceVolume[face0] + faceVolume[face1]);
1005 dA = faceVolume[face1] - faceVolume[face0];
1006 dfcc[0] = faceCentroid[face1][spaceId1] - faceCentroid[face0][spaceId1];
1007 dfcc[1] = faceCentroid[face1][spaceId2] - faceCentroid[face0][spaceId2];
1010 faceCentroid[face0][spaceId] + F1B3 * (F1 + Fvol * faceVolume[face1]) * cellLength[spaceId];
1014 * (faceVolume[face0] * faceCentroid[face0][spaceId1]
1015 + F1B2 * (faceVolume[face0] * dfcc[0] + dA * faceCentroid[face0][spaceId1]) + F1B3 * dA * dfcc[0]);
1019 * (faceVolume[face0] * faceCentroid[face0][spaceId2]
1020 + F1B2 * (faceVolume[face0] * dfcc[1] + dA * faceCentroid[face0][spaceId2]) + F1B3 * dA * dfcc[1]);
1022 for(
MInt i = 0; i < nDim; i++) {
1027 for(
MInt i = 0; i < nDim; i++) {
1040 for(
MInt i = 0; i < nDim; i++) {
1045 for(
MInt i = 0; i < nDim; i++) {
1051 * (faceVolume[0] * faceCentroid[0][0] - faceVolume[1] * faceCentroid[1][0]
1052 + faceVolume[2] * faceCentroid[2][1] - faceVolume[3] * faceCentroid[3][1]
1053 + faceVolume[4] * faceCentroid[4][2] - faceVolume[5] * faceCentroid[5][2]
1062 cerr <<
"this case has not yet been implemented" << endl;
1063 cerr <<
"cell " <<
cellId << endl;
1064 cerr <<
"coordinates"
1065 <<
" " << a_coordinate(cellId, 0) <<
" " << a_coordinate(cellId, 1) <<
" " << a_coordinate(cellId, 2)
1067 cerr <<
"number of cut faces: " << noCutFaces << endl;
1073 multimap<MFloat, MInt> sortedCPs;
1078 MFloat acentre[nDim] = {F0, F0, F0};
1080 for(
MInt i = 0; i < nDim; i++) {
1084 for(
MInt i = 0; i < nDim; i++) {
1087 for(
MInt i = 0; i < nDim; i++) {
1091 MFloat maxC = fabs(normal[0]);
1092 for(
MInt i = 1; i < nDim; i++) {
1093 if(fabs(normal[i]) > maxC) {
1094 maxC = fabs(normal[i]);
1098 spaceId1 = (spaceId + 1) % nDim;
1099 spaceId2 = (spaceId1 + 1) % nDim;
1100 vec_a[spaceId1] = F1;
1101 vec_a[spaceId2] = F1;
1102 vec_a[spaceId] = -(vec_a[spaceId1] * normal[spaceId1] + vec_a[spaceId2] * normal[spaceId2]) / normal[spaceId];
1104 for(
MInt i = 0; i < nDim; i++) {
1107 vec_b[spaceId] = normal[spaceId1] * vec_a[spaceId2] - normal[spaceId2] * vec_a[spaceId1];
1108 vec_b[spaceId1] = normal[spaceId2] * vec_a[spaceId] - normal[spaceId] * vec_a[spaceId2];
1109 vec_b[spaceId2] = normal[spaceId] * vec_a[spaceId1] - normal[spaceId1] * vec_a[spaceId];
1110 vecsum = sqrt(
POW2(vec_b[0]) +
POW2(vec_b[1]) +
POW2(vec_b[2]));
1111 for(
MInt i = 0; i < nDim; i++) {
1115 for(
MInt i = 0; i < nDim; i++) {
1120 for(
MInt i = 0; i < nDim; i++) {
1121 dx += vec_a[i] * (pCoords[i] - acentre[i]);
1122 dy += vec_b[i] * (pCoords[i] - acentre[i]);
1124 sortedCPs.insert(make_pair(atan2(dy, dx), cp));
1126 ASSERT((
signed)sortedCPs.size() == cutCell.
noCutPoints,
"");
1130 for(
auto& sortedCP : sortedCPs) {
1137 for(
MInt face = 0; face < m_noDirs; face++) {
1141 for(
MInt p = 0;
p < 4;
p++) {
1148 if(cutCell.
cutEdges[cp] == edgeCode2[face][p]) {
1164 if(improvedCentroidComputation) {
1172 MFloat bscentroid[nDim] = {F0, F0, F0};
1180 for(
MInt i = 0; i < nDim; i++) {
1184 crossProduct(vec0, vec1, vec2);
1186 for(
MInt i = 0; i < nDim; i++) {
1187 area +=
POW2(vec0[i]);
1189 area = F1B2 * sqrt(area);
1190 for(
MInt i = 0; i < nDim; i++) {
1197 for(
MInt i = 0; i < nDim; i++) {
1198 bscentroid[i] /=
mMax(1e-14, bsarea);
1200 for(
MInt i = 0; i < nDim; i++) {
1205 MFloat acentroid[nDim] = {F0, F0, F0};
1208 for(
MInt i = 0; i < nDim; i++) {
1214 for(
MInt i = 0; i < nDim; i++) {
1219 for(
MInt i = 0; i < nDim; i++) {
1220 acentroid[i] /= acnt;
1223 MFloat vcentroid[nDim] = {F0, F0, F0};
1229 for(
MInt i = 0; i < nDim; i++) {
1230 vec0[i] = a_coordinate(cellId, i) + signStencil[p0][i] * cellHalfLength[0] - acentroid[i];
1234 for(
MInt i = 0; i < nDim; i++) {
1235 vec0[i] = cutCell.
cutPoints[cp0][i] - acentroid[i];
1243 for(
MInt i = 0; i < nDim; i++) {
1244 vec1[i] = a_coordinate(cellId, i) + signStencil[p1][i] * cellHalfLength[0] - acentroid[i];
1248 for(
MInt i = 0; i < nDim; i++) {
1249 vec1[i] = cutCell.
cutPoints[cp1][i] - acentroid[i];
1253 for(
MInt i = 0; i < nDim; i++) {
1254 vec2[i] = a_coordinate(cellId, i) + signStencil[p2][i] * cellHalfLength[0] - acentroid[i];
1258 for(
MInt i = 0; i < nDim; i++) {
1259 vec2[i] = cutCell.
cutPoints[cp2][i] - acentroid[i];
1262 crossProduct(vec_tmp, vec0, vec1);
1264 for(
MInt i = 0; i < nDim; i++) {
1265 tvol += F1B6 * vec_tmp[i] * vec2[i];
1268 for(
MInt i = 0; i < nDim; i++) {
1270 fabs(tvol) * (F3B4 * (acentroid[i] + F1B3 * (vec0[i] + vec1[i] + vec2[i])) + F1B4 * acentroid[i]);
1274 for(
MInt i = 0; i < nDim; i++) {
1275 vcentroid[i] /= vol2;
1276 vcentroid[i] -= a_coordinate(cellId, i);
1279 for(
MInt i = 0; i < nDim; i++) {
1293 const MInt tCutGroup) {
1295 const MUint noCutCells = cutCellData.size();
1299 static constexpr MInt cornersSOLVERtoMC[8] = {1, 2, 0, 3, 5, 6, 4, 7};
1300 static constexpr MInt edgesMCtoSOLVER[12] = {0, 2, 1, 3, 4, 6, 5, 7, 10, 8, 9, 11};
1301 static constexpr MInt facesMCtoSOLVER[6] = {0, 2, 1, 3, 4, 5};
1302 static constexpr MFloat signStencil[8][3] = {{-F1, -F1, -F1}, {F1, -F1, -F1}, {-F1, F1, -F1}, {F1, F1, -F1},
1303 {-F1, -F1, F1}, {F1, -F1, F1}, {-F1, F1, F1}, {F1, F1, F1}};
1304 static constexpr MInt edgeCornerCode[12][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}, {4, 6}, {5, 7},
1305 {4, 5}, {6, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
1306 static constexpr MInt faceEdgeCode[6][4] = {{0, 10, 4, 8}, {1, 9, 5, 11}, {2, 8, 6, 9},
1307 {3, 11, 7, 10}, {2, 1, 3, 0}, {6, 4, 7, 5}};
1308 static constexpr MInt edgeFaceCode[12][2] = {{0, 4}, {1, 4}, {2, 4}, {3, 4}, {0, 5}, {1, 5},
1309 {2, 5}, {3, 5}, {0, 2}, {1, 2}, {0, 3}, {1, 3}};
1310 static constexpr MInt faceCornerCode[6][4] = {{0, 2, 6, 4}, {3, 1, 5, 7}, {1, 0, 4, 5},
1311 {2, 3, 7, 6}, {0, 1, 3, 2}, {5, 4, 6, 7}};
1312 const MInt maxNoSets = 6;
1313 ASSERT(m_noLevelSetsUsedForMb < maxNoSets,
"");
1315 MIntScratchSpace presentCases(m_noLevelSetsUsedForMb, 15, AT_,
"presentCases");
1316 for(
MInt s = 0; s < m_noLevelSetsUsedForMb; s++) {
1317 for(
MInt i = 0; i < 15; i++) {
1318 presentCases(s, i) = 0;
1322 const MInt maxNoFaces = 100;
1323 const MInt maxNoVertices = 300;
1324 const MInt maxNoEdges = 300;
1326 std::vector<polyVertex> vertices[maxNoSets];
1327 std::vector<polyFace> faces[maxNoSets];
1328 std::vector<polyEdge> edges[maxNoSets];
1329 std::vector<polyCutCell> cutCells;
1331 stack<MInt> faceStack;
1333 std::vector<polyMultiFace> multiFaces;
1334 MIntScratchSpace edgeCutCellPointer(maxNoEdges, AT_,
"edgeCutCellPointer");
1335 std::list<MInt> pureBodyEdges;
1337 MIntScratchSpace multiFaceConnection(maxNoFaces, AT_,
"multiFaceConnection");
1339 MIntScratchSpace outcode_MC_set(m_noLevelSetsUsedForMb, AT_,
"outcode_MC_set");
1340 MIntScratchSpace noCutPointsFromSet(m_noLevelSetsUsedForMb, AT_,
"cutSetPointers");
1341 MIntScratchSpace cutSetPointers(m_noLevelSetsUsedForMb, AT_,
"cutSetPointers");
1343 std::vector<CsgPolygon> csg_polygons;
1344 std::vector<CsgVertex> csg_vertices;
1345 std::vector<Csg> csg;
1346 std::vector<CsgPolygon> result;
1347 std::vector<polyVertex> vertices_result;
1348 MIntScratchSpace vertices_renamed(
mMax(2, m_noLevelSetsUsedForMb), maxNoVertices, AT_,
"vertices_renamed");
1349 std::list<std::pair<MInt, MInt>> openEdges;
1352 std::vector<MInt> faceVertices;
1355 "normalDotProduct");
1358 MBoolScratchSpace isPartOfBodySurface(maxNoVertices, AT_,
"isPartOfBodySurface");
1364 MIntScratchSpace surfaceIdentificatorCounters(m_noDirs + m_noEmbeddedBodies, AT_,
"surfaceIdentificatorCounters");
1365 MIntScratchSpace cellSurfaceMapping_backup(m_noDirs, AT_,
"cellSurfaceMapping_backup");
1368 NEW_SUB_TIMER_STATIC(tCutFace_0,
"cutFaceNew_0", tCutGroup);
1369 NEW_SUB_TIMER_STATIC(tCutFace_1a,
"cutFaceNew_1a", tCutGroup);
1370 NEW_SUB_TIMER_STATIC(tCutFace_1,
"cutFaceNew_1", tCutGroup);
1371 NEW_SUB_TIMER_STATIC(tCutFace_2,
"cutFaceNew_2", tCutGroup);
1372 NEW_SUB_TIMER_STATIC(tCutFace_3,
"cutFaceNew_3", tCutGroup);
1373 NEW_SUB_TIMER_STATIC(tCutFace_3b,
"cutFaceNew_3b", tCutGroup);
1374 NEW_SUB_TIMER_STATIC(tCutFace_4,
"cutFaceNew_4", tCutGroup);
1375 NEW_SUB_TIMER_STATIC(tCutFace_5a,
"cutFaceNew_5a", tCutGroup);
1376 NEW_SUB_TIMER_STATIC(tCutFace_5b,
"cutFaceNew_5b", tCutGroup);
1377 NEW_SUB_TIMER_STATIC(tCutFace_6,
"cutFaceNew_6", tCutGroup);
1378 NEW_SUB_TIMER_STATIC(tCutFace_7,
"cutFaceNew_7", tCutGroup);
1381 const MInt debugDomainId = 5;
1382 const MInt debugTimeStep = 0;
1383 const MInt debugCellId = 21866;
1392 const MFloat difBound1 = 10 * m_eps;
1393 const MFloat difBound2 = 100000000 * m_eps;
1396 for(
MUint cutc = 0; cutc < noCutCells; cutc++) {
1398 RECORD_TIMER_START(tCutFace_0);
1403 const MInt cellId = cutCellData[cutc].cellId;
1406 const MFloat cellLength0 = a_cellLengthAtCell(cellId);
1407 const MFloat cellHalfLength = F1B2 * cellLength0;
1409 for(
MBool& externalFace : cutCellData[cutc].externalFaces) {
1410 externalFace =
false;
1422 for(
MInt s = 0; s < m_noLevelSetsUsedForMb; s++) {
1423 vertices[s].clear();
1426 cutSetPointers[s] = -1;
1427 noCutPointsFromSet[s] = 0;
1430 const MBool isGapCell = cutCellData[cutc].isGapCell;
1434 if(m_complexBoundary && m_noLevelSetsUsedForMb > 1 && (!isGapCell)) {
1436 endSet = m_noLevelSetsUsedForMb;
1442 MBool isCompletelyOutside =
false;
1443 MBool hasAmbiguousFaces =
false;
1444 for(
MInt s = startSet; s < endSet; s++) {
1448 unsigned char outcode_MC = 0;
1449 for(
MInt c = 0; c < m_noCorners; c++) {
1451 MBool currentOutcode = (!cutCellData[cutc].cornerIsInsideGeometry[s][c]);
1452 if(currentOutcode) {
1453 outcode_MC = outcode_MC | (1 << cornersSOLVERtoMC[c]);
1456 outcode_MC_set[s] = outcode_MC;
1457 if(outcode_MC == 0) {
1458 isCompletelyOutside =
true;
1462 for(
MInt cutPoint = 0; cutPoint < cutCellData[cutc].noCutPoints; cutPoint++) {
1464 if(m_complexBoundary && m_noLevelSetsUsedForMb > 1 && (!isGapCell)) {
1466 set = m_bodyToSetTable[cutCellData[cutc].cutBodyIds[cutPoint]];
1469 noCutPointsFromSet[set]++;
1471 for(
MInt set = startSet; set < endSet; set++) {
1472 if(noCutPointsFromSet[set] && !isCompletelyOutside) {
1473 cutSetPointers[set] = noCutSets++;
1476 MInt maxCutsPerFace = 0;
1477 MInt cutsPerFace[m_noCorners] = {0, 0, 0, 0, 0, 0};
1478 for(
MInt cutPoint = 0; cutPoint < cutCellData[cutc].noCutPoints; cutPoint++) {
1479 MInt edge = cutCellData[cutc].cutEdges[cutPoint];
1480 cutsPerFace[edgeFaceCode[edge][0]]++;
1481 cutsPerFace[edgeFaceCode[edge][1]]++;
1483 for(
MInt k = 0; k < CC::noFaces; k++) {
1484 maxCutsPerFace =
mMax(maxCutsPerFace, cutsPerFace[k]);
1487 RECORD_TIMER_STOP(tCutFace_0);
1490 for(
MInt set = startSet; set < endSet; set++) {
1491 MInt setIndex = cutSetPointers[set];
1495 MInt outcode_MC = outcode_MC_set[set];
1497 MInt currentCase = cases[outcode_MC][0];
1498 if(!caseStatesLs[currentCase][0]) {
1499 cerr <<
"Error: Case not implemented: " << currentCase << endl;
1500 mTerm(1, AT_,
"Case not implemented, see error file for deatails.");
1503 if(!caseStatesLs[currentCase][1]) {
1504 hasAmbiguousFaces =
true;
1510 MBool simpleMarchingCubesSucceeds =
false;
1512 if(!hasAmbiguousFaces && noCutSets == 1 && maxCutsPerFace == 2 && !isCompletelyOutside
1513 && m_noLevelSetsUsedForMb == 1) {
1515 RECORD_TIMER_START(tCutFace_1a);
1517 simpleMarchingCubesSucceeds = computeCutFaceSimple(cutCellData[cutc]);
1518 RECORD_TIMER_STOP(tCutFace_1a);
1521 if(simpleMarchingCubesSucceeds) {
1526 for(
MInt i = 0; i < m_noLevelSetsUsedForMb; i++) {
1532 for(
MInt set = startSet; set < endSet; set++) {
1533 MInt setIndex = cutSetPointers[set];
1539 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
1540 cerr <<
"Cutting-Cell with " << set <<
" Index " << setIndex <<
" bodyId "
1541 << cutCellData[cutc].associatedBodyIds[set] << endl;
1545 RECORD_TIMER_START(tCutFace_1);
1548 MInt bodyId = cutCellData[cutc].associatedBodyIds[set];
1549 if(bodyId < 0)
continue;
1554 MInt cutPoints[12] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
1555 MInt cutPointToVertexMap[12] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
1556 MInt noCutPoints = 0;
1558 for(
MInt cutPoint = 0; cutPoint < cutCellData[cutc].noCutPoints; cutPoint++) {
1559 if(m_complexBoundary && m_noLevelSetsUsedForMb > 1 && (!isGapCell)) {
1561 if(m_bodyToSetTable[cutCellData[cutc].cutBodyIds[cutPoint]] != set) {
1567 cutPoints[cutCellData[cutc].cutEdges[cutPoint]] = cutPoint;
1569 vertices[setIndex].emplace_back(cutCellData[cutc].cutPoints[cutPoint], cutPoint, 1);
1570 cutPointToVertexMap[cutPoint] = vertices[setIndex].size() - 1;
1572 MInt cornerToVertexMap[8] = {-1, -1, -1, -1, -1, -1, -1, -1};
1573 for(
MInt c = 0; c < m_noCorners; c++) {
1575 if(!cutCellData[cutc].cornerIsInsideGeometry[set][c]) {
1576 std::array<MFloat, nDim> tmp_coords{};
1577 for(
MInt i = 0; i < nDim; i++) {
1578 tmp_coords[i] = a_coordinate(cellId, i) + signStencil[c][i] * cellHalfLength;
1580 vertices[setIndex].emplace_back(tmp_coords, c, 0);
1581 cornerToVertexMap[c] = vertices[setIndex].size() - 1;
1591 MInt outcode_MC = outcode_MC_set[set];
1594 MInt currentCase = cases[outcode_MC][0];
1595 presentCases(set, currentCase)++;
1596 MInt currentSubCase = cases[outcode_MC][1];
1598 if(!caseStatesLs[currentCase][0]) {
1599 cerr <<
"Error:Case not implemented: " << currentCase << endl;
1600 mTerm(1, AT_,
"Case not implemented, see error file for deatails.");
1604 m_caseCheckList[outcode_MC] = 1;
1608 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
1609 cerr <<
"cases: " << outcode_MC <<
" case " << currentCase <<
" subcase " << currentSubCase <<
" ambiguous "
1610 << caseStatesLs[currentCase][1] << endl;
1615 if(!caseStatesLs[currentCase][1]) {
1617 MInt faceMax = noAmbiguousFaces[currentCase];
1618 const MInt* testPointer;
1619 switch(currentCase) {
1621 testPointer = test3[currentSubCase];
1624 testPointer = test6[currentSubCase];
1627 testPointer = test7[currentSubCase];
1630 testPointer = test10[currentSubCase];
1633 testPointer = test12[currentSubCase];
1636 testPointer = test13[currentSubCase];
1639 mTerm(1, AT_,
"this should not happen!");
1642 for(
MInt i = 0; i < faceMax; i++) {
1643 MInt face = testPointer[i];
1644 MBool revert =
false;
1653 face = facesMCtoSOLVER[face];
1655 MBool testResult = cutCellData[cutc].faceCentroidIsInsideGeometry[set][face];
1657 if(revert) testResult = !testResult;
1660 subConfig +=
IPOW2(i);
1667 const MInt* tilingPointer =
nullptr;
1668 MBool addVertex =
false;
1669 noTriangles[set] = noTriangles_simpleCases[currentCase];
1670 switch(currentCase) {
1675 tilingPointer = tiling1Ls[currentSubCase];
1678 tilingPointer = tiling2Ls[currentSubCase];
1681 tilingPointer = tiling4Ls[currentSubCase];
1684 tilingPointer = tiling5Ls[currentSubCase];
1687 tilingPointer = tiling8Ls[currentSubCase];
1690 tilingPointer = tiling9Ls[currentSubCase];
1693 tilingPointer = tiling11Ls[currentSubCase];
1696 tilingPointer = tiling14Ls[currentSubCase];
1702 tilingPointer = tiling3_1Ls[currentSubCase];
1703 noTriangles[set] = 2;
1706 tilingPointer = tiling3_2Ls[currentSubCase];
1707 noTriangles[set] = 4;
1710 mTerm(1, AT_,
"This should not happen - we have case 3 with subConfig > 1... exiting");
1718 tilingPointer = tiling6_1_1Ls[currentSubCase];
1719 noTriangles[set] = 3;
1722 tilingPointer = tiling6_2Ls[currentSubCase];
1723 noTriangles[set] = 5;
1726 mTerm(1, AT_,
"This should not happen - we have case 6 with subConfig > 1... exiting");
1734 tilingPointer = tiling7_1Ls[currentSubCase];
1735 noTriangles[set] = 3;
1738 tilingPointer = tiling7_2Ls[currentSubCase][0];
1739 noTriangles[set] = 5;
1742 tilingPointer = tiling7_2Ls[currentSubCase][1];
1743 noTriangles[set] = 5;
1746 tilingPointer = tiling7_3Ls[currentSubCase][0];
1747 noTriangles[set] = 9;
1751 tilingPointer = tiling7_2Ls[currentSubCase][2];
1752 noTriangles[set] = 5;
1755 tilingPointer = tiling7_3Ls[currentSubCase][1];
1756 noTriangles[set] = 9;
1760 tilingPointer = tiling7_3Ls[currentSubCase][2];
1761 noTriangles[set] = 9;
1765 tilingPointer = tiling7_4_1Ls[currentSubCase];
1766 noTriangles[set] = 5;
1769 mTerm(1, AT_,
"This should not happen - we have case 7 with subConfig > 7... exiting");
1777 tilingPointer = tiling10_1_1Ls[currentSubCase];
1778 noTriangles[set] = 4;
1781 tilingPointer = tiling10_2Ls[currentSubCase];
1782 noTriangles[set] = 8;
1786 tilingPointer = tiling10_2_Ls[currentSubCase];
1787 noTriangles[set] = 8;
1791 tilingPointer = tiling10_1_1_Ls[currentSubCase];
1792 noTriangles[set] = 4;
1795 mTerm(1, AT_,
"This should not happen - we have case 10 with subConfig > 3... exiting");
1803 tilingPointer = tiling12_1_1Ls[currentSubCase];
1804 noTriangles[set] = 4;
1807 tilingPointer = tiling12_2Ls[currentSubCase];
1808 noTriangles[set] = 8;
1812 tilingPointer = tiling12_2_Ls[currentSubCase];
1813 noTriangles[set] = 8;
1817 tilingPointer = tiling12_1_1_Ls[currentSubCase];
1818 noTriangles[set] = 4;
1821 mTerm(1, AT_,
"This should not happen - we have case 12 with subConfig > 3... exiting");
1827 MInt config_identificator = 0;
1828 subConfig = subconfig13[subConfig];
1829 config_identificator = 0;
1832 tilingPointer = tiling13_1Ls[currentSubCase];
1833 noTriangles[set] = 4;
1842 config_identificator = subConfig - 1;
1843 tilingPointer = tiling13_2Ls[currentSubCase][config_identificator];
1844 noTriangles[set] = 6;
1859 config_identificator = subConfig - 7;
1860 tilingPointer = tiling13_3Ls[currentSubCase][config_identificator];
1861 noTriangles[set] = 10;
1869 config_identificator = subConfig - 19;
1870 tilingPointer = tiling13_4Ls[currentSubCase][config_identificator];
1871 noTriangles[set] = 12;
1879 config_identificator = subConfig - 23;
1880 tilingPointer = tiling13_5_1Ls[currentSubCase][config_identificator];
1881 noTriangles[set] = 6;
1896 config_identificator = subConfig - 27;
1897 tilingPointer = tiling13_3_Ls[currentSubCase][config_identificator];
1898 noTriangles[set] = 10;
1908 config_identificator = subConfig - 39;
1909 tilingPointer = tiling13_2_Ls[currentSubCase][config_identificator];
1910 noTriangles[set] = 6;
1914 tilingPointer = tiling13_1_Ls[currentSubCase];
1915 noTriangles[set] = 4;
1919 mTerm(1, AT_,
"impossible MC case 13 subConfig, how could this happen? exiting...");
1925 mTerm(1, AT_,
"invalid MC case, how could this happen? exiting...");
1928 RECORD_TIMER_STOP(tCutFace_1);
1931 RECORD_TIMER_START(tCutFace_2);
1933 MInt additionalVertexId = -1;
1935 additionalVertexId = vertices[setIndex].size();
1936 vertices[setIndex].emplace_back(-1, 2);
1937 addPoint(&vertices[setIndex], cutPointToVertexMap, noCutPoints,
1938 vertices[setIndex][additionalVertexId].coordinates);
1941 ASSERT(noCutPoints >= nDim,
"");
1942 ASSERT(noCutPoints == caseCutPoints[currentCase], to_string(grid().tree().globalId(cellId)));
1944 for(
MInt t = 0; t < noTriangles[set]; t++) {
1947 MInt currentFace = faces[setIndex].size();
1948 faces[setIndex].emplace_back(-1, 1, bodyId);
1951 for(
MInt pt = 0; pt < 3; pt++) {
1952 MInt cutEdge = tilingPointer[t * 3 + pt];
1954 p[pt] = additionalVertexId;
1956 p[pt] = cutPointToVertexMap[cutPoints[edgesMCtoSOLVER[cutEdge]]];
1958 ASSERT(p[pt] >= 0, to_string(cutEdge) +
" " + to_string(isGapCell) +
" " + to_string(cellId) +
" "
1959 + to_string(grid().tree().globalId(cellId)) +
" " + to_string(edgesMCtoSOLVER[cutEdge])
1960 +
" " + to_string(cutPoints[edgesMCtoSOLVER[cutEdge]]) +
" "
1961 + to_string(cutPointToVertexMap[cutPoints[edgesMCtoSOLVER[cutEdge]]]) +
" "
1962 + to_string(cutCellData[cutc].noCutPoints));
1964 for(
MInt te = 0; te < 3; te++) {
1966 ASSERT(point > -1,
"");
1967 vertices[setIndex][point].surfaceIdentificators.insert(bodyId + m_noDirs);
1968 MInt nextPoint = p[(te + 1) % 3];
1969 MBool edgeExists =
false;
1970 for(
MInt e = 0; (unsigned)e < vertices[setIndex][p[te]].edges.size(); e++) {
1971 MInt edge = vertices[setIndex][point].edges[e];
1972 if(edges[setIndex][edge].vertices[1] == nextPoint) {
1974 faces[setIndex][currentFace].edges.emplace_back(edge, 1);
1975 const MInt vertexId = edges[setIndex][edge].vertices[0];
1976 faces[setIndex][currentFace].vertices.push_back(vertexId);
1977 edges[setIndex][edge].face[0] = currentFace;
1979 }
else if(edges[setIndex][edge].vertices[0] == nextPoint) {
1981 faces[setIndex][currentFace].edges.emplace_back(edge, -1);
1982 const MInt vertexId = edges[setIndex][edge].vertices[1];
1983 faces[setIndex][currentFace].vertices.push_back(vertexId);
1985 edges[setIndex][edge].face[1] = currentFace;
1991 MInt newEdge = edges[setIndex].size();
1993 if(vertices[setIndex][point].
pointType == 2 || vertices[setIndex][nextPoint].
pointType == 2) {
1998 MInt p0 = vertices[setIndex][point].pointId;
1999 MInt p1 = vertices[setIndex][nextPoint].pointId;
2001 MInt edge0 = cutCellData[cutc].cutEdges[p0];
2003 MInt edge1 = cutCellData[cutc].cutEdges[p1];
2004 if(edgeFaceCode[edge0][0] == edgeFaceCode[edge1][0] || edgeFaceCode[edge0][0] == edgeFaceCode[edge1][1])
2005 edgeId = edgeFaceCode[edge0][0];
2006 else if(edgeFaceCode[edge0][1] == edgeFaceCode[edge1][0]
2007 || edgeFaceCode[edge0][1] == edgeFaceCode[edge1][1])
2008 edgeId = edgeFaceCode[edge0][1];
2010 edges[setIndex].emplace_back(point, nextPoint, edgeId, edgeType);
2011 vertices[setIndex][point].edges.push_back(newEdge);
2012 vertices[setIndex][nextPoint].edges.push_back(newEdge);
2013 faces[setIndex][currentFace].edges.emplace_back(newEdge, 1);
2014 const MInt vertexId = edges[setIndex][newEdge].vertices[0];
2015 faces[setIndex][currentFace].vertices.push_back(vertexId);
2017 edges[setIndex][newEdge].face[0] = currentFace;
2021 computeNormal(vertices[setIndex][p[0]].coordinates, vertices[setIndex][p[1]].coordinates,
2022 vertices[setIndex][p[2]].coordinates, faces[setIndex][currentFace].normal,
2023 faces[setIndex][currentFace].w);
2026 for(
auto& f : faces[setIndex]) {
2027 ASSERT(f.vertices.size() == 3,
"");
2031 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2032 cerr <<
"Body surfaces: " << endl;
2033 cerr <<
"Set: " << set <<
" noTriangles " << noTriangles[set] <<
" noFaces " << faces[setIndex].size()
2034 <<
" noEdges " << edges[setIndex].size() <<
" noVertices " << vertices[setIndex].size() << endl;
2039 RECORD_TIMER_STOP(tCutFace_2);
2041 RECORD_TIMER_START(tCutFace_3);
2048 for(
MInt e = 0; e < m_noEdges; e++) {
2050 MBool p0Fluid = !cutCellData[cutc].cornerIsInsideGeometry[set][edgeCornerCode[e][0]];
2052 MBool p1Fluid = !cutCellData[cutc].cornerIsInsideGeometry[set][edgeCornerCode[e][1]];
2053 if(p0Fluid && p1Fluid) {
2054 MInt v0 = cornerToVertexMap[edgeCornerCode[e][0]];
2055 MInt v1 = cornerToVertexMap[edgeCornerCode[e][1]];
2056 edges[setIndex].emplace_back(v0, v1, e, 0);
2057 vertices[setIndex][v0].edges.push_back(edges[setIndex].size() - 1);
2058 vertices[setIndex][v1].edges.push_back(edges[setIndex].size() - 1);
2059 }
else if(p0Fluid) {
2060 MInt v0 = cornerToVertexMap[edgeCornerCode[e][0]];
2061 MInt v1 = cutPointToVertexMap[cutPoints[e]];
2062 edges[setIndex].emplace_back(v0, v1, e, 1);
2063 vertices[setIndex][v0].edges.push_back(edges[setIndex].size() - 1);
2064 vertices[setIndex][v1].edges.push_back(edges[setIndex].size() - 1);
2065 }
else if(p1Fluid) {
2066 MInt v0 = cornerToVertexMap[edgeCornerCode[e][1]];
2067 MInt v1 = cutPointToVertexMap[cutPoints[e]];
2068 edges[setIndex].emplace_back(v0, v1, e, 1);
2069 vertices[setIndex][v0].edges.push_back(edges[setIndex].size() - 1);
2070 vertices[setIndex][v1].edges.push_back(edges[setIndex].size() - 1);
2076 for(
MInt cartFace = 0; cartFace < m_noDirs; cartFace++) {
2078 MBool edgeTouched[4] = {
false,
false,
false,
false};
2079 for(
MInt e = 0; e < 4; e++) {
2080 if(edgeTouched[e])
continue;
2081 MInt cartEdge = faceEdgeCode[cartFace][e];
2082 MInt startVertex = cornerToVertexMap[faceCornerCode[cartFace][e]];
2084 if(startVertex == -1)
continue;
2087 MInt currentFace = faces[setIndex].size();
2088 faces[setIndex].emplace_back(cartFace, 0, -1);
2089 for(
MInt i = 0; i < nDim; i++) {
2090 faces[setIndex][currentFace].normal[i] = F0;
2092 MInt normalDir = cartFace / 2;
2094 if(cartFace % 2 == 0) {
2097 faces[setIndex][currentFace].w = -vertices[setIndex][startVertex].coordinates[normalDir] * normalSign;
2098 faces[setIndex][currentFace].normal[normalDir] = normalSign;
2102 MInt currentVertex = startVertex;
2104 vertices[setIndex][currentVertex].surfaceIdentificators.insert(cartFace);
2105 edgeTouched[currentE] =
true;
2108 MBool edgeFound =
false;
2109 for(
MInt ve = 0; (unsigned)ve < vertices[setIndex][currentVertex].edges.size(); ve++) {
2110 MInt edge = vertices[setIndex][currentVertex].edges[ve];
2111 if(edges[setIndex][edge].edgeType > 1)
2113 if(edges[setIndex][edge].edgeId != cartEdge)
continue;
2116 MInt edgeDirection = 1;
2117 if(edges[setIndex][edge].vertices[1] == currentVertex) {
2120 faces[setIndex][currentFace].edges.emplace_back(edge, edgeDirection);
2121 MInt vertexId = edges[setIndex][edge].vertices[0];
2122 if(edgeDirection == -1) vertexId = edges[setIndex][edge].vertices[1];
2123 faces[setIndex][currentFace].vertices.push_back(vertexId);
2125 if(edges[setIndex][edge].face[0] > -1)
2126 edges[setIndex][edge].face[1] = currentFace;
2128 edges[setIndex][edge].face[0] = currentFace;
2129 if(edgeDirection == 1)
2130 currentVertex = edges[setIndex][edge].vertices[1];
2132 currentVertex = edges[setIndex][edge].vertices[0];
2135 ASSERT(edgeFound,
"");
2138 if(vertices[setIndex][currentVertex].
pointType == 0) {
2139 currentE = (currentE + 1) % 4;
2140 cartEdge = faceEdgeCode[cartFace][currentE];
2141 }
else if(vertices[setIndex][currentVertex].
pointType == 1) {
2143 for(
MInt ve = 0; (unsigned)ve < vertices[setIndex][currentVertex].edges.size(); ve++) {
2144 MInt edge = vertices[setIndex][currentVertex].edges[ve];
2145 if(edges[setIndex][edge].edgeType != 2)
2147 if(edges[setIndex][edge].edgeId != cartFace)
continue;
2150 MInt edgeDirection = 1;
2151 if(edges[setIndex][edge].vertices[1] == currentVertex) {
2154 faces[setIndex][currentFace].edges.emplace_back(edge, edgeDirection);
2155 MInt vertexId = edges[setIndex][edge].vertices[0];
2156 if(edgeDirection == -1) vertexId = edges[setIndex][edge].vertices[1];
2157 faces[setIndex][currentFace].vertices.push_back(vertexId);
2159 edges[setIndex][edge].face[1] = currentFace;
2160 if(edgeDirection == 1)
2161 currentVertex = edges[setIndex][edge].vertices[1];
2163 currentVertex = edges[setIndex][edge].vertices[0];
2165 cartEdge = cutCellData[cutc].cutEdges[vertices[setIndex][currentVertex].pointId];
2166 while(faceEdgeCode[cartFace][currentE] != cartEdge) {
2167 currentE = (currentE + 1) % 4;
2171 ASSERT(edgeFound,
"");
2173 }
while(currentVertex != startVertex);
2179 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2180 cerr <<
"Including cartesian faces: " << endl;
2181 cerr <<
"Set: " << set <<
" SetIndex:" << setIndex <<
" faces : " << faces[setIndex].size()
2182 <<
" edges : " << edges[setIndex].size() <<
" certices : " << vertices[setIndex].size() << endl;
2183 for(
MUint i = 0; i < vertices[setIndex].size(); i++) {
2184 cerr <<
" vertices " << setprecision(19) << vertices[setIndex][i].coordinates[0] <<
" "
2185 << vertices[setIndex][i].coordinates[1] <<
" " << vertices[setIndex][i].coordinates[2] << endl;
2191 RECORD_TIMER_STOP(tCutFace_3);
2194 RECORD_TIMER_START(tCutFace_3b);
2196 MInt noIndividualCuts = 0;
2197 for(
MInt set = startSet; set < endSet; set++) {
2198 if(cutSetPointers[set] > -1) {
2204 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2205 cerr <<
"noIndividualCuts: " << noIndividualCuts << endl;
2209 if(noIndividualCuts == 0) {
2210 cutCellData[cutc].volume = F0;
2211 RECORD_TIMER_STOP(tCutFace_3b);
2214 ASSERT(noIndividualCuts > 0,
"");
2215 ASSERT(noCutSets > 0,
" ");
2216 ASSERT(noCutSets == noIndividualCuts,
"");
2220 MBool error =
false;
2221 const MInt startSetIndex = 0;
2222 MInt referenceSet = startSet;
2223 for(
MInt set = startSet; set < endSet; set++) {
2224 if(cutSetPointers[set] > -1) {
2237 for(
MInt set = startSet; set < endSet; set++) {
2238 MInt setIndex = cutSetPointers[set];
2239 if(setIndex < 0)
continue;
2241 for(
MInt faceId = 0; faceId < (signed)faces[setIndex].size(); faceId++) {
2242 ASSERT(faces[setIndex][faceId].edges.size() == faces[setIndex][faceId].vertices.size(),
"");
2243 for(
MInt edgeId = 0; edgeId < (signed)faces[setIndex][faceId].edges.size(); edgeId++) {
2244 const MInt edge = faces[setIndex][faceId].edges[edgeId].first;
2245 const MInt newVertexId = faces[setIndex][faceId].vertices[edgeId];
2246 const MInt direction = faces[setIndex][faceId].edges[edgeId].second;
2247 MInt vertexId = edges[setIndex][edge].vertices[0];
2248 if(direction == -1) vertexId = edges[setIndex][edge].vertices[1];
2249 ASSERT(newVertexId == vertexId,
"");
2258 MInt noInitialFaces = 0;
2264 if(noIndividualCuts > 1) {
2266 ASSERT(m_complexBoundary && m_noLevelSetsUsedForMb > 1,
"");
2267 ASSERT(m_bodyFaceJoinMode == 4 || m_bodyFaceJoinMode == 1,
"ERROR: using unrecommended bodyFaceJoinMode");
2271 for(
MInt set = referenceSet + 1; set < endSet; set++) {
2272 MInt setIndex = cutSetPointers[set];
2273 if(setIndex < 0)
continue;
2276 csg_polygons.clear();
2280 for(
MInt face = faces[startSetIndex].size() - 1; face > -1; face--) {
2281 csg_vertices.clear();
2282 polyFace* faceP = &faces[startSetIndex][face];
2284 for(
MInt e = 0; (unsigned)e < faceP->vertices.size(); e++) {
2286 csg_vertices.emplace_back(
CsgVector(vertices[startSetIndex][vertexId].coordinates), vertexId,
2292 csg_polygons.emplace_back(csg_vertices, startSetIndex, faceP->
faceId, faceP->
faceType, faceP->
bodyId,
2295 csg_polygons.emplace_back(csg_vertices, startSetIndex, faceP->
faceId, faceP->
faceType, faceP->
bodyId);
2299 csg.emplace_back(csg_polygons);
2300 csg_polygons.clear();
2303 for(
MInt face = faces[setIndex].size() - 1; face > -1; face--) {
2304 csg_vertices.clear();
2305 polyFace* faceP = &faces[setIndex][face];
2307 for(
MInt e = 0; (unsigned)e < faceP->vertices.size(); e++) {
2309 csg_vertices.emplace_back(
CsgVector(vertices[setIndex][vertexId].coordinates), vertexId, setIndex);
2314 csg_polygons.emplace_back(csg_vertices, setIndex, faceP->
faceId, faceP->
faceType, faceP->
bodyId, plane);
2316 csg_polygons.emplace_back(csg_vertices, setIndex, faceP->
faceId, faceP->
faceType, faceP->
bodyId);
2319 csg.emplace_back(csg_polygons);
2322 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2323 cerr <<
"Combine: reference " << referenceSet <<
" with " << set <<
"Index " << startSetIndex <<
"face-size "
2324 << faces[startSetIndex].size() <<
"Index " << setIndex <<
"face-size " << faces[setIndex].size() << endl;
2332 result = csg[0].intersect(csg[1]);
2333 vertices_result.clear();
2337 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2338 cerr <<
"result-size: " << result.size() <<
"vertices_result " << vertices_result.size() << endl;
2340 for(
MInt poligon = 0; (unsigned)poligon < result.size(); poligon++) {
2341 cerr <<
"resultId " << poligon <<
" body " << result[poligon].bodyId <<
" type " << result[poligon].faceType
2342 <<
" setIndex " << setIndex << endl;
2343 for(
MInt vertex = 0; (unsigned)vertex < result[poligon].vertices.size(); vertex++) {
2344 cerr <<
"vertex " << vertex <<
"coordinate: " << setprecision(21)
2345 << result[poligon].vertices[vertex].pos.xx[0] <<
" " << result[poligon].vertices[vertex].pos.xx[1]
2346 <<
" " << result[poligon].vertices[vertex].pos.xx[2] <<
" "
2347 << result[poligon].vertices[vertex].vertexId <<
" " << result[poligon].vertices[vertex].setIndex
2350 cerr <<
"faceId " << result[poligon].faceId << endl;
2362 for(
MInt i = 0; i <
mMax(2, m_noLevelSetsUsedForMb); i++) {
2363 for(
MInt j = 0; j < maxNoVertices; j++) {
2364 vertices_renamed(i, j) = -1;
2369 MInt noVertices = 0;
2375 for(
MInt p = 0; (unsigned)p < result.size(); p++) {
2376 ASSERT(result[p].vertices.size() <= (
unsigned)maxNoVertices,
"");
2377 ASSERT((
signed)result[p].vertices.size() >= 3,
"");
2380 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2381 MInt vertexId = result[p].vertices[v].vertexId;
2382 MInt vertexSetIndex = result[p].vertices[v].setIndex;
2383 if(vertexSetIndex < 0)
continue;
2385 if(vertexId < 0)
continue;
2389 if(vertices_renamed(vertexSetIndex, vertexId) == -1) {
2391 MBool vertexFound =
false;
2392 MInt vertexIndex = -1;
2395 for(
MInt i = 0; (unsigned)i < vertices_result.size(); i++) {
2397 coord_diff += pow(result[p].vertices[v].pos.xx[0] - vertices_result[i].coordinates[0], 2);
2398 coord_diff += pow(result[p].vertices[v].pos.xx[1] - vertices_result[i].coordinates[1], 2);
2399 coord_diff += pow(result[p].vertices[v].pos.xx[2] - vertices_result[i].coordinates[2], 2);
2401 if(m_multiCutCell) coord_diff = pow(coord_diff, 0.5);
2404 if(coord_diff < difBound1) {
2406 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2407 cerr <<
"Small diff " << p <<
" " << vertexId <<
" " << vertices_result[i].coordinates[0] <<
" "
2408 << vertices_result[i].coordinates[1] <<
" " << vertices_result[i].coordinates[2] <<
" eps "
2409 << difBound1 <<
" dif " << coord_diff <<
" " << result[p].vertices[v].pos.xx[0] <<
" "
2410 << result[p].vertices[v].pos.xx[1] <<
" " << result[p].vertices[v].pos.xx[2] <<
" " << endl;
2420 vertices_renamed(vertexSetIndex, vertexId) = vertexIndex;
2423 vertices_renamed(vertexSetIndex, vertexId) = noVertices;
2424 vertices_result.emplace_back(vertices[vertexSetIndex][vertexId].pointId,
2425 vertices[vertexSetIndex][vertexId].
pointType);
2427 ASSERT(vertices[vertexSetIndex][vertexId].
pointType == 0
2428 || vertices[vertexSetIndex][vertexId].
pointType == 1
2429 || vertices[vertexSetIndex][vertexId].
pointType == 2,
2431 vertices_result[noVertices].coordinates[0] = vertices[vertexSetIndex][vertexId].coordinates[0];
2432 vertices_result[noVertices].coordinates[1] = vertices[vertexSetIndex][vertexId].coordinates[1];
2433 vertices_result[noVertices].coordinates[2] = vertices[vertexSetIndex][vertexId].coordinates[2];
2438 result[p].vertices[v].vertexId = vertices_renamed(vertexSetIndex, vertexId);
2439 result[p].vertices[v].setIndex = -1;
2446 for(
MInt p = 0; (unsigned)p < result.size(); p++) {
2447 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2448 MInt vertexId = result[p].vertices[v].vertexId;
2449 if(vertexId == -1) {
2451 MBool vertexFound =
false;
2452 MInt vertexIndex = -1;
2453 for(
MInt i = 0; (unsigned)i < vertices_result.size(); i++) {
2455 coord_diff += pow(result[p].vertices[v].pos.xx[0] - vertices_result[i].coordinates[0], 2);
2456 coord_diff += pow(result[p].vertices[v].pos.xx[1] - vertices_result[i].coordinates[1], 2);
2457 coord_diff += pow(result[p].vertices[v].pos.xx[2] - vertices_result[i].coordinates[2], 2);
2459 if(m_multiCutCell) coord_diff = pow(coord_diff, 0.5);
2462 MFloat difBound = difBound1;
2463 if(m_multiCutCell) {
2464 difBound = difBound2;
2467 if(coord_diff < difBound) {
2469 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2470 cerr <<
"2 Small diff " << p <<
" " << vertexId <<
" " << vertices_result[i].coordinates[0] <<
" "
2471 << vertices_result[i].coordinates[1] <<
" " << vertices_result[i].coordinates[2] <<
" eps "
2472 << coord_diff <<
" " << result[p].vertices[v].pos.xx[0] <<
" "
2473 << result[p].vertices[v].pos.xx[1] <<
" " << result[p].vertices[v].pos.xx[2] <<
" " << endl;
2483 result[p].vertices[v].vertexId = vertexIndex;
2484 result[p].vertices[v].setIndex = -1;
2488 vertices_result.emplace_back(-1, 3);
2489 vertices_result[noVertices].coordinates[0] = result[p].vertices[v].pos.xx[0];
2490 vertices_result[noVertices].coordinates[1] = result[p].vertices[v].pos.xx[1];
2491 vertices_result[noVertices].coordinates[2] = result[p].vertices[v].pos.xx[2];
2492 result[p].vertices[v].vertexId = noVertices;
2493 result[p].vertices[v].setIndex = -1;
2501 vertices[startSetIndex].swap(vertices_result);
2502 vertices_result.clear();
2503 edges[startSetIndex].clear();
2504 faces[startSetIndex].clear();
2507 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2508 cerr <<
"VertexId-check: " << endl;
2510 for(
MInt poligon = 0; (unsigned)poligon < result.size(); poligon++) {
2511 cerr <<
"resultId " << poligon <<
" body " << result[poligon].bodyId <<
" type " << result[poligon].faceType
2512 <<
" setIndex " << setIndex << endl;
2513 for(
MInt vertex = 0; (unsigned)vertex < result[poligon].vertices.size(); vertex++) {
2514 cerr <<
"vertex " << vertex <<
"coordinate: " << setprecision(21)
2515 << result[poligon].vertices[vertex].pos.xx[0] <<
" " << result[poligon].vertices[vertex].pos.xx[1]
2516 <<
" " << result[poligon].vertices[vertex].pos.xx[2] <<
" "
2517 << result[poligon].vertices[vertex].vertexId <<
" " << result[poligon].vertices[vertex].setIndex
2520 cerr <<
"faceId " << result[poligon].faceId << endl;
2523 for(
MInt poligon = 0; (unsigned)poligon < result.size(); poligon++) {
2524 for(
MInt vertex = 0; (unsigned)vertex < result[poligon].vertices.size(); vertex++) {
2525 MInt vertexId = result[poligon].vertices[vertex].vertexId;
2526 for(
MInt i = 0; i < nDim; i++) {
2527 if(fabs(result[poligon].vertices[vertex].pos.xx[i] - vertices[startSetIndex][vertexId].coordinates[i])
2529 cerr <<
" VertexId-missmatch " << cellId <<
" " << vertexId <<
" " << vertex <<
" " << i
2530 << setprecision(21) << result[poligon].vertices[vertex].pos.xx[i] <<
" "
2531 << vertices[startSetIndex][vertexId].coordinates[i] << endl;
2549 for(
MInt p = 0; (unsigned)p < result.size(); p++) {
2551 MBool polygonValid =
true;
2552 MInt validVertices = 0;
2553 faceMapping[p] = -1;
2555 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2556 const MInt j = (v + 1) % result[p].vertices.
size();
2557 const MInt vertexId = result[p].vertices[v].vertexId;
2558 const MInt vertexIdNext = result[p].vertices[j].vertexId;
2559 if(vertexId == vertexIdNext) {
2560 vertexValid[v] =
false;
2562 vertexValid[v] =
true;
2566 if(vertexValid[v]) {
2567 MInt surfaceIndicator = -1;
2568 if(result[p].faceType == 0) {
2569 surfaceIndicator = result[p].faceId;
2571 surfaceIndicator = result[p].bodyId + m_noDirs;
2573 vertices[startSetIndex][vertexId].surfaceIdentificators.insert(surfaceIndicator);
2576 if(validVertices < 3) {
2590 polygonValid =
false;
2593 if(!polygonValid)
continue;
2595 faces[startSetIndex].emplace_back(result[p].faceId, result[p].faceType, result[p].bodyId);
2596 faces[startSetIndex][noFaces].normal[0] = result[p].plane.normal.xx[0];
2597 faces[startSetIndex][noFaces].normal[1] = result[p].plane.normal.xx[1];
2598 faces[startSetIndex][noFaces].normal[2] = result[p].plane.normal.xx[2];
2600 !(std::isnan(result[p].plane.normal.xx[0] + result[p].plane.normal.xx[1] + result[p].plane.normal.xx[2])),
2602 faces[startSetIndex][noFaces].tmpSetIndex = result[p].setIndex;
2603 faces[startSetIndex][noFaces].w = -result[p].plane.w;
2604 faceMapping[p] = noFaces;
2605 faces[startSetIndex][noFaces].edges.clear();
2606 faces[startSetIndex][noFaces].vertices.clear();
2609 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2610 const MInt vertexId = result[p].vertices[v].vertexId;
2611 MBool vertexUsed =
false;
2612 for(
MInt vertice : faces[startSetIndex][noFaces].vertices) {
2613 if(vertexId == vertice) {
2620 faces[startSetIndex][noFaces].vertices.push_back(vertexId);
2621 vertices[startSetIndex][vertexId].faceIds.push_back(noFaces);
2624 ASSERT((
signed)faces[startSetIndex][noFaces].vertices.size() >= 3,
"");
2630 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2631 cerr <<
"Current-Face-Count " << noFaces << endl;
2638#if defined CutCell_DEBUG || !defined NDEBUG
2639 MBool onlySingleEdges =
true;
2641 for(
MInt p = 0; (unsigned)p < result.size(); p++) {
2642 MInt faceId = faceMapping[p];
2643 if(faceId < 0)
continue;
2646 for(
MInt v = 0; (unsigned)v < result[p].vertices.size(); v++) {
2647 MInt j = (v + 1) % result[p].vertices.
size();
2648 MInt vertexId = result[p].vertices[v].vertexId;
2649 MInt vertexIdNext = result[p].vertices[j].vertexId;
2650 if(vertexId == vertexIdNext) {
2655 MBool edgeFound =
false;
2657 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
2658 MInt v0 = edges[startSetIndex][e].vertices[0];
2659 MInt v1 = edges[startSetIndex][e].vertices[1];
2660 if(vertexId == v0 && vertexIdNext == v1) {
2666 if(vertexId == v1 && vertexIdNext == v0) {
2675 edges[startSetIndex].emplace_back(vertexId, vertexIdNext, -1, -1);
2676 edges[startSetIndex][noEdges].face[0] = faceId;
2679 vertices[startSetIndex][vertexId].edges.push_back(noEdges);
2680 vertices[startSetIndex][vertexIdNext].edges.push_back(noEdges);
2684 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2685 cerr <<
"Adding Edge " << noEdges <<
" " << vertexId <<
" " << vertexIdNext <<
" " << faceId << endl;
2690 }
else if(edges[startSetIndex][edgeId].face[1] < 0) {
2691 edges[startSetIndex][edgeId].face[1] = faceId;
2694 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2695 cerr <<
"Adding 2nd face to egde " << edgeId <<
" " << vertexId <<
" " << vertexIdNext <<
" " << faceId
2697 cerr <<
"edge vertices are: " << edges[startSetIndex][edgeId].vertices[0] <<
" "
2698 << edges[startSetIndex][edgeId].vertices[1] << endl;
2704#if defined CutCell_DEBUG || !defined NDEBUG
2705 onlySingleEdges =
false;
2707 edges[startSetIndex].emplace_back(vertexId, vertexIdNext, -1, -1);
2708 edges[startSetIndex][noEdges].face[0] = faceId;
2711 vertices[startSetIndex][vertexId].edges.push_back(noEdges);
2712 vertices[startSetIndex][vertexIdNext].edges.push_back(noEdges);
2715 cerr << cellId <<
"Has a duplicate edge!" << endl;
2716 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2717 cerr <<
"Adding Duplicate Edge " << noEdges <<
" " << vertexId <<
" " << vertexIdNext <<
" " << faceId
2722 faces[startSetIndex][faceId].edges.emplace_back(edgeId, direction);
2727 for(
MInt faceId = 0; faceId < static_cast<MInt>(faces[startSetIndex].size()); faceId++) {
2728 ASSERT(faces[startSetIndex][faceId].edges.size() == faces[startSetIndex][faceId].vertices.size(),
"");
2729 for(
MInt edgeId = 0; edgeId < static_cast<MInt>(faces[startSetIndex][faceId].edges.size()); edgeId++) {
2730 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
2731 const MInt newVertexId = faces[startSetIndex][faceId].vertices[edgeId];
2732 const MInt direction = faces[startSetIndex][faceId].edges[edgeId].second;
2733 MInt vertexId = edges[startSetIndex][edge].vertices[0];
2734 if(direction == -1) {
2735 vertexId = edges[startSetIndex][edge].vertices[1];
2737 ASSERT(newVertexId == vertexId,
"");
2749 for(
MInt e = 0; e < static_cast<MInt>(edges[startSetIndex].size()); e++) {
2750 ASSERT(edges[startSetIndex][e].face[0] > -1,
"");
2751 if(edges[startSetIndex][e].face[1] == -1) {
2752 MInt face = edges[startSetIndex][e].face[0];
2754 for(
MInt fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
2755 if(faces[startSetIndex][face].edges[fe].first == e) {
2756 dir = faces[startSetIndex][face].edges[fe].second;
2761 openEdges.emplace_back(e, -1);
2763 openEdges.emplace_back(e, 1);
2768 vector<std::pair<MInt, MInt>> openEdgeList;
2770 openEdgeId.
fill(-1);
2771 for(
auto& openEdge : openEdges) {
2772 openEdgeId(openEdge.first) =
static_cast<MInt>(openEdgeList.size());
2773 openEdgeList.push_back(openEdge);
2777 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
2778 cerr << openEdges.
size() <<
" initially open edges! " << endl;
2779 cerr <<
"Count: " << faces[startSetIndex].size() <<
" faces " << edges[startSetIndex].size() <<
" edges "
2780 << vertices[startSetIndex].size() <<
" vertices " << endl;
2782 for(
auto it2 = openEdges.begin(); it2 != openEdges.end(); it2++) {
2783 const MInt edgdeId = (*it2).first;
2784 cerr <<
" openEdge with Id " << (*it2).first <<
" v1 " << edges[startSetIndex][edgdeId].vertices[0] <<
" "
2785 << edges[startSetIndex][edgdeId].vertices[1] << endl;
2799#if defined CutCell_DEBUG || !defined NDEBUG
2800 for(
MInt i = 0; i < (signed)openEdgeList.size(); i++) {
2801 const MInt e = openEdgeList[i].first;
2802 const MInt dir = openEdgeList[i].second;
2803 const MInt id0 = (dir == 1) ? 0 : 1;
2804 const MInt id1 = (dir == 1) ? 1 : 0;
2805 const MInt v0 = edges[startSetIndex][e].vertices[id0];
2806 const MInt v1 = edges[startSetIndex][e].vertices[id1];
2808 ASSERT(edges[startSetIndex][e].face[0] > -1,
"");
2811 for(
MInt j = 0; j < (signed)vertices[startSetIndex][v0].edges.size(); j++) {
2812 const MInt e2 = vertices[startSetIndex][v0].edges[j];
2814 if(e2 == e)
continue;
2816 const MInt face = edges[startSetIndex][e2].face[0];
2818 for(
MInt fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
2819 if(faces[startSetIndex][face].edges[fe].first == e2) {
2820 dir2 = faces[startSetIndex][face].edges[fe].second;
2823 const MInt id00 = (dir2 == 1) ? 0 : 1;
2824 const MInt id11 = (dir2 == 1) ? 1 : 0;
2825 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
2826 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
2828 ASSERT(edges[startSetIndex][e2].face[0] > -1,
"");
2829 if(onlySingleEdges) {
2830 ASSERT(v0 != v00 || v1 != v11,
"");
2831 ASSERT(v0 != v11 || v1 != v00,
"");
2840 MBool somethingChanged = !openEdgeList.empty();
2841 while(somethingChanged) {
2842 somethingChanged =
false;
2844 MBool alreadyResolved =
false;
2846 for(
MUint k1 = 0; k1 < openEdgeList.size(); k1++) {
2848 const MInt e1 = openEdgeList[k1].first;
2849 if(openEdgeId(e1) < 0)
continue;
2850 const MInt dir1 = openEdgeList[k1].second;
2851 const MInt rdir1 = (dir1 == 1) ? -1 : 1;
2852 const MInt id0 = (dir1 == 1) ? 0 : 1;
2853 const MInt id1 = (dir1 == 1) ? 1 : 0;
2854 const MInt v0 = edges[startSetIndex][e1].vertices[id0];
2855 const MInt v1 = edges[startSetIndex][e1].vertices[id1];
2856 const MInt face = edges[startSetIndex][e1].face[0];
2860 for(fe = 0; (unsigned)fe < faces[startSetIndex][face].edges.size(); fe++) {
2861 if(faces[startSetIndex][face].edges[fe].first == e1) {
2865 ASSERT(fe < (
signed)faces[startSetIndex][face].edges.size(),
"");
2870 for(fv = 0; (unsigned)fv < faces[startSetIndex][face].vertices.size(); fv++) {
2871 if(faces[startSetIndex][face].vertices[fv] == v0 || faces[startSetIndex][face].vertices[fv] == v1) {
2872 if(fv + 1 < (
signed)faces[startSetIndex][face].vertices.size()
2873 && (faces[startSetIndex][face].vertices[fv + 1] == v0
2874 || faces[startSetIndex][face].vertices[fv + 1] == v1)) {
2883 for(
MInt i = 0; i < nDim; i++) {
2884 a[i] = vertices[startSetIndex][v1].coordinates[i] - vertices[startSetIndex][v0].coordinates[i];
2885 a_abs +=
a[i] *
a[i];
2887 a_abs = sqrt(a_abs);
2890 for(
MInt k2 = 0; (unsigned)k2 < vertices[startSetIndex][v0].edges.size(); k2++) {
2892 const MInt e2 = vertices[startSetIndex][v0].edges[k2];
2893 if(e2 == e1)
continue;
2894 if(openEdgeId(e2) < 0)
continue;
2895 ASSERT(openEdgeList[openEdgeId(e2)].first == e2,
"");
2897 const MInt dir2 = openEdgeList[openEdgeId(e2)].second;
2898 const MInt id00 = (dir2 == 1) ? 0 : 1;
2899 const MInt id11 = (dir2 == 1) ? 1 : 0;
2900 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
2901 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
2903 if(v11 != v0)
continue;
2908 for(
MInt i = 0; i < nDim; i++) {
2909 b[i] = vertices[startSetIndex][v11].coordinates[i] - vertices[startSetIndex][v00].coordinates[i];
2910 b_abs +=
b[i] *
b[i];
2911 dotProduct +=
a[i] *
b[i];
2913 b_abs = sqrt(b_abs);
2915 dotProduct = fabs(F1 + (dotProduct / (a_abs * b_abs)));
2918 if(dotProduct < 1e-10 && b_abs < a_abs) {
2932 MBool foundMatchingEdge =
false;
2933 for(
MUint k3 = 0; k3 < openEdgeList.size(); k3++) {
2934 e3 = openEdgeList[k3].first;
2935 if(openEdgeId(e3) < 0)
continue;
2936 dir3 = openEdgeList[k3].second;
2937 id000 = (dir3 == 1) ? 0 : 1;
2938 id111 = (dir3 == 1) ? 1 : 0;
2939 v000 = edges[startSetIndex][e3].vertices[id000];
2940 v111 = edges[startSetIndex][e3].vertices[id111];
2941 if(v111 == v00 && v000 == v1) {
2942 foundMatchingEdge =
true;
2945 if(v000 == v00 && v111 == v1) {
2946 foundMatchingEdge =
true;
2971 edges[startSetIndex][e1].vertices[id0] = v00;
2974 edges[startSetIndex][e2].face[1] = face;
2977 MInt otherDir = (edges[startSetIndex][e2].vertices[id0] == v0) ? rdir1 : dir1;
2978 auto pos = faces[startSetIndex][face].edges.begin() + fe;
2979 pos += (dir1 == 1) ? 0 : 1;
2980 faces[startSetIndex][face].edges.insert(pos, make_pair(e2, otherDir));
2983 if(!m_multiCutCell) {
2985 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
2986 if(vertices[startSetIndex][v0].edges[edg] == e1) {
2987 vertices[startSetIndex][v0].edges[edg] = e2;
2992 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
2993 if(vertices[startSetIndex][v0].edges[edg] == e1) {
2994 auto epos = vertices[startSetIndex][v0].edges.begin() + edg;
2995 vertices[startSetIndex][v0].edges.erase(epos);
3000 if(m_multiCutCell) vertices[startSetIndex][v00].edges.push_back(e1);
3003 if(foundMatchingEdge) {
3004 MInt otherDir2 = (edges[startSetIndex][e2].vertices[id0] == v0) ? rdir1 : dir1;
3007 for(fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
3008 if(faces[startSetIndex][face].edges[fe].first == e1) {
3009 faces[startSetIndex][face].edges[fe].first = e3;
3010 faces[startSetIndex][face].edges[fe].second = otherDir2;
3015 edges[startSetIndex][e3].face[1] = face;
3018 edges[startSetIndex][e1].face[0] = -1;
3019 edges[startSetIndex][e1].face[1] = -1;
3022 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3023 cerr <<
"Replacing edge " << e1 <<
" with " << e3 <<
" and adding vertex " << v00 <<
" to face "
3025 cerr <<
"Intermediate edge is " << e2 <<
" with v " << v1 <<
" v2 " << v0 << endl;
3029 if(!m_multiCutCell) {
3031 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v00].edges.size(); edg++) {
3032 if(vertices[startSetIndex][v00].edges[edg] == e1) {
3033 vertices[startSetIndex][v00].edges[edg] = e3;
3036 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3037 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3038 vertices[startSetIndex][v1].edges[edg] = e3;
3043 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v00].edges.size(); edg++) {
3044 if(vertices[startSetIndex][v00].edges[edg] == e1) {
3045 auto epos = vertices[startSetIndex][v00].edges.begin() + edg;
3046 vertices[startSetIndex][v00].edges.erase(epos);
3049 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3050 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3051 auto epos = vertices[startSetIndex][v1].edges.begin() + edg;
3052 vertices[startSetIndex][v1].edges.erase(epos);
3058 openEdgeId(e1) = -1;
3059 openEdgeId(e3) = -1;
3063 auto vpos = faces[startSetIndex][face].vertices.
begin() + fv;
3064 if(vpos > faces[startSetIndex][face].vertices.end()) {
3065 faces[startSetIndex][face].vertices.push_back(v00);
3067 faces[startSetIndex][face].vertices.insert(vpos, v00);
3071 openEdgeId(e2) = -1;
3073 somethingChanged =
true;
3074 alreadyResolved =
true;
3079 if(alreadyResolved)
continue;
3082 for(
MInt k2 = 0; (unsigned)k2 < vertices[startSetIndex][v1].edges.size(); k2++) {
3083 const MInt e2 = vertices[startSetIndex][v1].edges[k2];
3084 if(e2 == e1)
continue;
3085 if(openEdgeId(e2) < 0)
continue;
3086 ASSERT(openEdgeList[openEdgeId(e2)].first == e2,
"");
3087 const MInt dir2 = openEdgeList[openEdgeId(e2)].second;
3088 const MInt id00 = (dir2 == 1) ? 0 : 1;
3089 const MInt id11 = (dir2 == 1) ? 1 : 0;
3090 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
3091 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
3093 if(v00 != v1)
continue;
3098 for(
MInt i = 0; i < nDim; i++) {
3099 b[i] = vertices[startSetIndex][v11].coordinates[i] - vertices[startSetIndex][v00].coordinates[i];
3100 b_abs +=
b[i] *
b[i];
3101 dotProduct +=
a[i] *
b[i];
3103 b_abs = sqrt(b_abs);
3105 dotProduct = fabs(F1 + (dotProduct / (a_abs * b_abs)));
3108 if(dotProduct < 1e-10 && b_abs < a_abs) {
3122 MBool foundMatchingEdge =
false;
3123 for(
auto& k3 : openEdgeList) {
3125 if(openEdgeId(e3) < 0)
continue;
3127 id000 = (dir3 == 1) ? 0 : 1;
3128 id111 = (dir3 == 1) ? 1 : 0;
3129 v000 = edges[startSetIndex][e3].vertices[id000];
3130 v111 = edges[startSetIndex][e3].vertices[id111];
3131 if(v000 == v0 && v111 == v11) {
3132 foundMatchingEdge =
true;
3134 }
else if(v000 == v11 && v111 == v0) {
3135 foundMatchingEdge =
true;
3159 MInt otherDir = (edges[startSetIndex][e2].vertices[id1] == v1) ? rdir1 : dir1;
3160 auto pos = faces[startSetIndex][face].edges.begin() + fe;
3161 pos += (dir1 == 1) ? 1 : 0;
3162 faces[startSetIndex][face].edges.insert(pos, make_pair(e2, otherDir));
3165 auto vpos = faces[startSetIndex][face].vertices.begin() + fv;
3166 if(vpos > faces[startSetIndex][face].vertices.end()) {
3167 faces[startSetIndex][face].vertices.push_back(v11);
3169 faces[startSetIndex][face].vertices.insert(vpos, v11);
3173 edges[startSetIndex][e2].face[1] = face;
3176 if(!m_multiCutCell) {
3177 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3178 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3179 vertices[startSetIndex][v1].edges[edg] = e2;
3184 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v1].edges.size(); edg++) {
3185 if(vertices[startSetIndex][v1].edges[edg] == e1) {
3186 auto epos = vertices[startSetIndex][v1].edges.begin() + edg;
3187 vertices[startSetIndex][v1].edges.erase(epos);
3192 if(m_multiCutCell) vertices[startSetIndex][v11].edges.push_back(e1);
3195 edges[startSetIndex][e1].vertices[id1] = v11;
3198 openEdgeId(e2) = -1;
3201 if(foundMatchingEdge) {
3205 const MInt otherDir2 = (edges[startSetIndex][e3].vertices[id1] == v1) ? rdir1 : dir1;
3208 for(fe = 0; fe < (signed)faces[startSetIndex][face].edges.size(); fe++) {
3209 if(faces[startSetIndex][face].edges[fe].first == e1) {
3210 faces[startSetIndex][face].edges[fe].first = e3;
3211 faces[startSetIndex][face].edges[fe].second = otherDir2;
3217 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3218 cerr <<
"Replacing edge " << e1 <<
" with " << e3 <<
" and adding vertex " << v00 <<
" to face "
3220 cerr <<
"Intermediate edge is " << e2 <<
" with v " << v11 <<
" v2 " << v1 << endl;
3225 edges[startSetIndex][e1].face[0] = -1;
3226 edges[startSetIndex][e1].face[1] = -1;
3229 if(!m_multiCutCell) {
3230 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
3231 if(vertices[startSetIndex][v0].edges[edg] == e1) {
3232 vertices[startSetIndex][v0].edges[edg] = e3;
3235 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v11].edges.size(); edg++) {
3236 if(vertices[startSetIndex][v11].edges[edg] == e1) {
3237 vertices[startSetIndex][v11].edges[edg] = e3;
3242 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v0].edges.size(); edg++) {
3243 if(vertices[startSetIndex][v0].edges[edg] == e1) {
3244 auto epos = vertices[startSetIndex][v0].edges.
begin() + edg;
3245 vertices[startSetIndex][v0].edges.erase(epos);
3248 for(
MInt edg = 0; (unsigned)edg < vertices[startSetIndex][v11].edges.size(); edg++) {
3249 if(vertices[startSetIndex][v11].edges[edg] == e1) {
3250 auto epos = vertices[startSetIndex][v11].edges.begin() + edg;
3251 vertices[startSetIndex][v11].edges.erase(epos);
3256 openEdgeId(e1) = -1;
3257 openEdgeId(e3) = -1;
3260 edges[startSetIndex][e3].face[1] = face;
3263 somethingChanged =
true;
3275 for(
MUint i = 0; i < openEdgeList.size(); i++) {
3276 const MInt e = openEdgeList[i].first;
3277 const MInt dir = openEdgeList[i].second;
3278 const MInt id0 = (dir == 1) ? 0 : 1;
3279 const MInt id1 = (dir == 1) ? 1 : 0;
3280 const MInt v0 = edges[startSetIndex][e].vertices[id0];
3281 const MInt v1 = edges[startSetIndex][e].vertices[id1];
3284 if(edges[startSetIndex][e].face[0] == -1 && edges[startSetIndex][e].face[1] == -1)
continue;
3287 for(
MInt j = 0; (unsigned)j < vertices[startSetIndex][v0].edges.size(); j++) {
3288 const MInt e2 = vertices[startSetIndex][v0].edges[j];
3289 if(e2 == e)
continue;
3290 const MInt dir2 = openEdgeList[openEdgeId(e2)].second;
3291 const MInt id00 = (dir2 == 1) ? 0 : 1;
3292 const MInt id11 = (dir2 == 1) ? 1 : 0;
3293 const MInt v00 = edges[startSetIndex][e2].vertices[id00];
3294 const MInt v11 = edges[startSetIndex][e2].vertices[id11];
3297 if(edges[startSetIndex][e2].face[0] == -1 && edges[startSetIndex][e2].face[1] == -1)
continue;
3298 if(onlySingleEdges) {
3299 ASSERT(v0 != v00 || v1 != v11,
"");
3300 ASSERT(v0 != v11 || v1 != v00,
"");
3307 for(
MInt faceId = 0; faceId < (signed)faces[startSetIndex].size(); faceId++) {
3308 for(
MInt edgeId = 0; edgeId < (signed)faces[startSetIndex][faceId].edges.size(); edgeId++) {
3309 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
3310 ASSERT(edges[startSetIndex][edge].face[0] > -1,
"");
3315 for(
MInt v = 0; v < (signed)vertices[startSetIndex].size(); v++) {
3316 for(
MInt edgeId = 0; (unsigned)edgeId < vertices[startSetIndex][v].edges.size(); edgeId++) {
3317 const MInt edge = vertices[startSetIndex][v].edges[edgeId];
3318 ASSERT(edges[startSetIndex][edge].face[0] > -1,
"");
3323 for(
MInt faceId = 0; faceId < (signed)faces[startSetIndex].size(); faceId++) {
3324 ASSERT(faces[startSetIndex][faceId].edges.size() == faces[startSetIndex][faceId].vertices.size(),
"");
3325 for(
MInt edgeId = 0; edgeId < (signed)faces[startSetIndex][faceId].edges.size(); edgeId++) {
3326 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
3327 const MInt newVertexId = faces[startSetIndex][faceId].vertices[edgeId];
3328 const MInt direction = faces[startSetIndex][faceId].edges[edgeId].second;
3329 MInt vertexId = edges[startSetIndex][edge].vertices[0];
3330 if(direction == -1) vertexId = edges[startSetIndex][edge].vertices[1];
3331 ASSERT(newVertexId == vertexId,
"");
3336 std::list<std::pair<MInt, MInt>> openEdgesOld(openEdges);
3338 for(
auto& it1 : openEdgesOld) {
3339 if(openEdgeId(it1.first) < 0)
continue;
3340 openEdges.push_back(it1);
3345 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3346 cerr << openEdges.size() <<
" open edge remaining: " << endl;
3347 cerr <<
"Count: " << faces[startSetIndex].size() <<
" faces " << edges[startSetIndex].size() <<
" edges "
3348 << vertices[startSetIndex].size() <<
" vertices " << endl;
3350 for(
auto it2 = openEdges.begin(); it2 != openEdges.end(); it2++) {
3351 const MInt edgeId = (*it2).first;
3352 cerr <<
" openEdge with Id " << edgeId <<
" v1 " << edges[startSetIndex][edgeId].vertices[0] <<
" "
3353 << edges[startSetIndex][edgeId].vertices[1] << endl;
3365 if(openEdges.size() == 1) {
3366 cerr <<
" Only 1 open edge remaining: This means trouble! " << endl;
3367 cerr <<
"Re-checking coordinate-distances! " << endl;
3368 for(
MInt i = 0; (unsigned)i < vertices[startSetIndex].size(); i++) {
3369 for(
MInt j = 0; (unsigned)j < vertices[startSetIndex].size(); j++) {
3370 if(i == j)
continue;
3373 pow(vertices[startSetIndex][i].coordinates[0] - vertices[startSetIndex][j].coordinates[0], 2);
3375 pow(vertices[startSetIndex][i].coordinates[1] - vertices[startSetIndex][j].coordinates[1], 2);
3377 pow(vertices[startSetIndex][i].coordinates[2] - vertices[startSetIndex][j].coordinates[2], 2);
3379 MFloat dif_old = coord_diff;
3380 coord_diff = pow(coord_diff, 0.5);
3382 if(coord_diff < cellLength0) {
3383 cerr <<
" Coord-dif " << coord_diff <<
" i " << i <<
" " << j << endl;
3384 cerr <<
" dif-old " << dif_old <<
" " << m_eps * 10 << endl;
3398 for(
auto it1 = openEdges.begin(); it1 != openEdges.end(); it1++) {
3400 const MInt e1 = (*it1).first;
3401 const MInt dir1 = (*it1).second;
3402 MInt v11 = edges[startSetIndex][e1].vertices[0];
3403 MInt v12 = edges[startSetIndex][e1].vertices[1];
3411 for(
MInt i = 0; i < nDim; i++) {
3412 a[i] = vertices[startSetIndex][v12].coordinates[i] - vertices[startSetIndex][v11].coordinates[i];
3413 a_abs +=
a[i] *
a[i];
3415 a_abs = sqrt(a_abs);
3417 for(
auto& openEdge : openEdges) {
3418 const MInt e2 = openEdge.first;
3419 const MInt dir2 = openEdge.second;
3422 dotProdMatrix(e1, e2) = 1000.0;
3423 dotProdMatrix(e2, e1) = 1000.0;
3426 MInt v21 = edges[startSetIndex][e2].vertices[0];
3427 MInt v22 = edges[startSetIndex][e2].vertices[1];
3439 for(
MInt i = 0; i < nDim; i++) {
3440 b[i] = vertices[startSetIndex][v22].coordinates[i] - vertices[startSetIndex][v21].coordinates[i];
3441 b_abs +=
b[i] *
b[i];
3442 dotProduct +=
a[i] *
b[i];
3444 b_abs = sqrt(b_abs);
3445 dotProduct = abs(F1 - abs(dotProduct / (a_abs * b_abs)));
3446 dotProdMatrix(e1, e2) = dotProduct;
3449 dotProdMatrix(e1, e2) = 1000.0;
3456 while(!openEdges.empty()) {
3457 faceVertices.clear();
3460 MInt eLast = openEdges.front().first;
3461 MInt dir = openEdges.front().second;
3465 for(
MInt i = 0; i < m_noDirs + m_noEmbeddedBodies; i++) {
3466 surfaceIdentificatorCounters[i] = 0;
3470 faces[startSetIndex].emplace_back(faceId, faceType, bodyId);
3472 faces[startSetIndex][noFaces].edges.emplace_back(eLast, dir);
3473 faces[startSetIndex][noFaces].isLine = 1;
3476 openEdges.pop_front();
3478 MInt startVertex = edges[startSetIndex][eLast].vertices[0];
3479 MInt vertex = edges[startSetIndex][eLast].vertices[1];
3482 MInt tmp = startVertex;
3483 startVertex = vertex;
3487 faces[startSetIndex][noFaces].vertices.push_back(startVertex);
3490 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3491 cerr <<
"Adding face line starting with edge " << eLast <<
" with vertex " << startVertex <<
" and "
3497 faceVertices.push_back(vertex);
3499 for(
MInt surfaceIdentificator : vertices[startSetIndex][vertex].surfaceIdentificators) {
3500 surfaceIdentificatorCounters[surfaceIdentificator]++;
3503 edges[startSetIndex][eLast].face[1] = noFaces;
3506 while(vertex != startVertex) {
3508 auto bestEdgeIt = openEdges.
end();
3509 MFloat minDotProduct = 1001.0;
3514 for(
auto it = openEdges.begin(); it != openEdges.end(); it++) {
3517 if(dotProdMatrix(eLast, e) < minDotProduct) {
3519 dir = (*bestEdgeIt).second;
3520 MInt vStart = edges[startSetIndex][e].vertices[0];
3521 MInt vEnd = edges[startSetIndex][e].vertices[1];
3527 if(vStart == vertex) {
3528 minDotProduct = dotProdMatrix(eLast, e);
3530 }
else if(vEnd == vertex) {
3531 minDotProduct = dotProdMatrix(eLast, e);
3537 if(bestEdgeIt == openEdges.end()) {
3538 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex],
globalTimeStep);
3539 cerr << grid().domainId() <<
"open edges for cell " << cellId <<
": ";
3540 cerr <<
"Open edges size: " << openEdges.size() << endl;
3541 for(
auto& openEdge : openEdges) {
3542 cerr << openEdge.first;
3544 mTerm(1, AT_,
"Open edge for cutCell!");
3548 e = (*bestEdgeIt).first;
3549 dir = (*bestEdgeIt).second;
3550 MInt vStart = edges[startSetIndex][e].vertices[0];
3551 MInt vEnd = edges[startSetIndex][e].vertices[1];
3557 if(vStart != vertex) {
3558 if(vEnd == vertex) {
3565 MInt vEndOld = vEnd;
3569 cerr << grid().domainId() <<
" missmatching edges for " << cellId <<
" "
3570 << grid().tree().globalId(cellId) <<
" vertices: " << vStart <<
" " << vEnd <<
" " << vertex
3578 faces[startSetIndex][noFaces].edges.emplace_back(e, dir);
3579 faces[startSetIndex][noFaces].vertices.push_back(vStart);
3582 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3583 cerr <<
"Adding edge " << e <<
" with vertex " << vStart <<
" and " << vEnd << endl;
3584 cerr <<
"The new face now has " << faces[startSetIndex][noFaces].vertices.size() <<
" vertices!" << endl;
3588 openEdges.erase(bestEdgeIt);
3590 faceVertices.push_back(vertex);
3591 for(
MInt surfaceIdentificator : vertices[startSetIndex][vertex].surfaceIdentificators)
3592 surfaceIdentificatorCounters[surfaceIdentificator]++;
3593 edges[startSetIndex][e].face[1] = noFaces;
3596 if(faceVertices.size() <= 2) {
3597 cerr << grid().domainId() <<
" face with only two vertices for " << cellId <<
" "
3598 << grid().tree().globalId(cellId) <<
" face : " << noFaces << endl;
3602 MInt noDoubleVertices = 0;
3603 ASSERT(vertices[startSetIndex].size() <= (
unsigned)maxNoVertices,
"");
3604 for(
MInt v = 0; (unsigned)v < vertices[startSetIndex].size(); v++) {
3605 vertexTouches[v] = 0;
3607 for(
MInt fv = 0; (unsigned)fv < faceVertices.size(); fv++) {
3608 vertexTouches[faceVertices[fv]]++;
3610 for(
MInt v = 0; (unsigned)v < vertices[startSetIndex].size(); v++) {
3611 if(vertexTouches[v] > 1) {
3616 if(noDoubleVertices != 0) {
3617 cerr <<
"double vertices for " << cellId <<
" on " << grid().domainId() << endl;
3623 MInt bestIdentificator = -1;
3625 for(
MInt i = 0; i < m_noDirs + m_noEmbeddedBodies; i++) {
3626 if(surfaceIdentificatorCounters[i] > maxHits) {
3627 maxHits = surfaceIdentificatorCounters[i];
3628 bestIdentificator = i;
3631 ASSERT(bestIdentificator > -1,
"");
3632 ASSERT(bestIdentificator < m_noDirs + m_noEmbeddedBodies,
"");
3633 faceType = (bestIdentificator < m_noDirs ? 0 : 1);
3635 faceId = bestIdentificator;
3637 bodyId = bestIdentificator - m_noDirs;
3639 faces[startSetIndex][noFaces].faceType = faceType;
3640 faces[startSetIndex][noFaces].faceId = faceId;
3641 faces[startSetIndex][noFaces].bodyId = bodyId;
3644 MBool partnerFound =
false;
3645 for(
MInt e2 = 0; (unsigned)e2 < faces[startSetIndex][noFaces].edges.size(); e2++) {
3646 MInt edge = faces[startSetIndex][noFaces].edges[e2].first;
3647 MInt otherFace = edges[startSetIndex][edge].face[0];
3648 MInt otherFaceType = faces[startSetIndex][otherFace].faceType;
3649 MInt otherFaceId = faces[startSetIndex][otherFace].faceId;
3650 MInt otherBodyId = faces[startSetIndex][otherFace].bodyId;
3651 if(otherFaceType == faceType && otherFaceId == faceId && otherBodyId == bodyId) {
3652 partnerFound =
true;
3653 faces[startSetIndex][noFaces].normal[0] = faces[startSetIndex][otherFace].normal[0];
3654 faces[startSetIndex][noFaces].normal[1] = faces[startSetIndex][otherFace].normal[1];
3655 faces[startSetIndex][noFaces].normal[2] = faces[startSetIndex][otherFace].normal[2];
3656 faces[startSetIndex][noFaces].w = faces[startSetIndex][otherFace].w;
3662 ASSERT(partnerFound,
"");
3672 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3673 cerr <<
"Final count: " << faces[startSetIndex].size() <<
" faces " << edges[startSetIndex].size() <<
" edges "
3674 << vertices[startSetIndex].size() <<
" vertices " << endl;
3680 for(
MInt faceId = 0; faceId < (signed)faces[startSetIndex].size(); faceId++) {
3681 ASSERT(faces[startSetIndex][faceId].edges.size() == faces[startSetIndex][faceId].vertices.size(),
"");
3682 for(
MInt edgeId = 0; edgeId < (signed)faces[startSetIndex][faceId].edges.size(); edgeId++) {
3683 const MInt edge = faces[startSetIndex][faceId].edges[edgeId].first;
3684 const MInt newVertexId = faces[startSetIndex][faceId].vertices[edgeId];
3685 const MInt direction = faces[startSetIndex][faceId].edges[edgeId].second;
3686 MInt vertexId = edges[startSetIndex][edge].vertices[0];
3687 if(direction == -1) vertexId = edges[startSetIndex][edge].vertices[1];
3688 ASSERT(newVertexId == vertexId,
"");
3692 RECORD_TIMER_STOP(tCutFace_3b);
3693 RECORD_TIMER_START(tCutFace_4);
3697 ASSERT(faceStack.empty(),
"");
3700 if(faces[startSetIndex].size() == 0) {
3701 cutCells.emplace_back(cellId, &a_coordinate(cellId, 0));
3705 for(
MInt faceCounter = 0; (unsigned)faceCounter < faces[startSetIndex].size(); faceCounter++) {
3706 if(faces[startSetIndex][faceCounter].cutCell > -1)
continue;
3709 if(faces[startSetIndex][faceCounter].isLine)
continue;
3712 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3713 cerr <<
"faceCounter: " << faceCounter << endl;
3717 faceStack.push(faceCounter);
3718 const MInt currentCutCell = cutCells.size();
3719 cutCells.emplace_back(cellId, &a_coordinate(cellId, 0));
3720 faces[startSetIndex][faceCounter].cutCell = currentCutCell;
3721 cutCells[currentCutCell].faces.push_back(faceCounter);
3723 while(!faceStack.empty()) {
3724 MInt currentFace = faceStack.top();
3727 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][currentFace].edges.size(); e++) {
3728 MInt edge = faces[startSetIndex][currentFace].edges[e].first;
3729 MInt otherFace = edges[startSetIndex][edge].face[0];
3730 if(otherFace == currentFace) {
3731 otherFace = edges[startSetIndex][edge].face[1];
3735 if(faces[startSetIndex][otherFace].cutCell == -1) {
3737 if(grid().domainId() == debugDomainId &&
globalTimeStep == debugTimeStep && cellId == debugCellId) {
3738 cerr <<
"faceCounter: " << faceCounter <<
" cutCell " << currentCutCell
3739 <<
" is adding neighbor: " << otherFace << endl;
3742 cutCells[currentCutCell].faces.push_back(otherFace);
3743 faces[startSetIndex][otherFace].cutCell = currentCutCell;
3744 faceStack.push(otherFace);
3749 ASSERT(faceStack.empty(),
"");
3751 RECORD_TIMER_STOP(tCutFace_4);
3754 RECORD_TIMER_START(tCutFace_5a);
3755 compVolumeIntegrals_pyraBased3(&cutCells, &faces[startSetIndex], &vertices[startSetIndex]);
3756 RECORD_TIMER_STOP(tCutFace_5a);
3759 if(
globalTimeStep == debugTimeStep && cellId == debugCellId && grid().domainId() == debugDomainId) {
3760 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], -73);
3765 cerr <<
"[" << -1 <<
"]"
3766 <<
": Warning: removal/suppression of doubly touched vertices seems to have failed!" << endl;
3767 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], -1);
3771 RECORD_TIMER_START(tCutFace_5b);
3782 if(m_bodyFaceJoinMode != 0) {
3785 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
3786 MInt face = edges[startSetIndex][e].face[0];
3788 edgeCutCellPointer[e] = -1;
3790 edgeCutCellPointer[e] = faces[startSetIndex][face].cutCell;
3792 for(
MInt f = 0; (unsigned)f < faces[startSetIndex].size(); f++) {
3793 multiFaceConnection[f] = -1;
3795 for(
MInt c = 0; (unsigned)c < cutCells.size(); c++) {
3796 switch(m_bodyFaceJoinMode) {
3801 list<MInt> bodyVertices;
3802 set<MInt> locBodySrfcs;
3804 isBodyEdge.
fill(-2);
3807 for(
auto& vert : vertices[startSetIndex]) {
3808 sort(vert.edges.begin(), vert.edges.end());
3809 auto last = unique(vert.edges.begin(), vert.edges.end());
3810 vert.edges.erase(last, vert.edges.end());
3813 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
3814 if(edgeCutCellPointer[e] != c)
continue;
3815 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 0
3816 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 0) {
3817 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceId
3818 == faces[startSetIndex][edges[startSetIndex][e].face[1]].faceId) {
3819 pureBodyEdges.push_back(e);
3821 }
else if(m_bodyFaceJoinMode == 1 || m_bodyFaceJoinMode == 4) {
3822 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 1
3823 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 1) {
3824 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId
3825 == faces[startSetIndex][edges[startSetIndex][e].face[1]].bodyId) {
3826 pureBodyEdges.push_back(e);
3827 isBodyEdge(e) = faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId;
3828 bodyVertices.push_back(edges[startSetIndex][e].vertices[0]);
3829 bodyVertices.push_back(edges[startSetIndex][e].vertices[1]);
3830 locBodySrfcs.insert(faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId);
3838 if(m_bodyFaceJoinMode == 4 && noIndividualCuts > 1) {
3839 bodyVertices.sort();
3840 bodyVertices.unique();
3841 const MInt noLocBodySrfcs = (signed)locBodySrfcs.
size();
3842 MInt concaveCnt = 0;
3843 for(
MInt vert : bodyVertices) {
3844 if(vertices[startSetIndex][vert].
pointType > 1) {
3845 MInt noOuterEdges = 0;
3846 MInt noBodyEdges = 0;
3847 MInt outerEdges[6] = {-1};
3848 MInt bodyEdges[6] = {-1};
3849 for(
MInt edge : vertices[startSetIndex][vert].edges) {
3850 if(isBodyEdge(edge) == -1) {
3851 outerEdges[noOuterEdges] = edge;
3854 if(isBodyEdge(edge) > -1) {
3855 bodyEdges[noBodyEdges] = edge;
3859 if(noOuterEdges == 2 && noBodyEdges > 0) {
3860 for(
MInt k = 0; k < noBodyEdges; k++) {
3861 MInt edge = bodyEdges[k];
3862 MInt otherVert = (edges[startSetIndex][edge].vertices[0] == vert) ? 1 : 0;
3863 otherVert = edges[startSetIndex][edge].vertices[otherVert];
3864 MInt face0 = edges[startSetIndex][edge].face[0];
3865 MInt face1 = edges[startSetIndex][edge].face[1];
3871 MInt noEdges0 = (signed)faces[startSetIndex][face0].edges.size();
3872 MInt noEdges1 = (signed)faces[startSetIndex][face1].edges.size();
3876 for(
MInt e = 0; e < noEdges0; e++) {
3877 if(faces[startSetIndex][face0].edges[e].first == edge) {
3879 dir0 = faces[startSetIndex][face0].edges[e].second;
3880 if(prevEdge < 0 || nextEdge < 0) {
3881 MInt vA = edges[startSetIndex][edge].vertices[dir0 == 1 ? 1 : 0];
3882 MInt otherEdge = (vA == vert) ? (edge0 + 1) % noEdges0 : (edge0 + noEdges0 - 1) % noEdges0;
3883 otherEdge = faces[startSetIndex][face0].edges[otherEdge].first;
3884 if(otherEdge == outerEdges[0] || otherEdge == outerEdges[1]) {
3885 MInt otherEdge2 = (otherEdge == outerEdges[0]) ? outerEdges[1] : outerEdges[0];
3886 nextEdge = (vA == vert) ? otherEdge : otherEdge2;
3887 prevEdge = (vA == vert) ? otherEdge2 : otherEdge;
3893 for(
MInt e = 0; e < noEdges1; e++) {
3894 if(faces[startSetIndex][face1].edges[e].first == edge) {
3896 dir1 = faces[startSetIndex][face1].edges[e].second;
3897 if(prevEdge < 0 || nextEdge < 0) {
3898 MInt vA = edges[startSetIndex][edge].vertices[dir1 == 1 ? 1 : 0];
3899 MInt otherEdge = (vA == vert) ? (edge1 + 1) % noEdges1 : (edge1 + noEdges1 - 1) % noEdges1;
3900 otherEdge = faces[startSetIndex][face1].edges[otherEdge].first;
3901 if(otherEdge == outerEdges[0] || otherEdge == outerEdges[1]) {
3902 MInt otherEdge2 = (otherEdge == outerEdges[0]) ? outerEdges[1] : outerEdges[0];
3903 nextEdge = (vA == vert) ? otherEdge : otherEdge2;
3904 prevEdge = (vA == vert) ? otherEdge2 : otherEdge;
3910 if(edge0 < 0 || edge1 < 0)
mTerm(1, AT_,
"ERROR: pure body edge not found.");
3911 if(dir0 == dir1)
mTerm(1, AT_,
"ERROR: pure body edge not found (2).");
3912 if(faces[startSetIndex][face0].edges[edge0].first != edge
3913 || faces[startSetIndex][face1].edges[edge1].first != edge)
3914 mTerm(1, AT_,
"ERROR: pure body edge not found (3).");
3915 ASSERT(outerEdges[0] > -1 && outerEdges[1] > -1,
"");
3917 if(prevEdge > -1 && nextEdge > -1) {
3918 ASSERT(side > -1,
"");
3919 MInt id_p0 = edges[startSetIndex][prevEdge].vertices[0] == vert ? 0 : 1;
3920 MInt id_p1 = edges[startSetIndex][nextEdge].vertices[0] == vert ? 0 : 1;
3921 if(edges[startSetIndex][prevEdge].vertices[id_p0] != vert)
3922 mTerm(1, AT_,
"ERROR: pure body edge not found (4).");
3923 if(edges[startSetIndex][nextEdge].vertices[id_p1] != vert)
3924 mTerm(1, AT_,
"ERROR: pure body edge not found (5).");
3925 MInt otherId[2] = {1, 0};
3926 MInt vA2 = edges[startSetIndex][prevEdge].vertices[otherId[id_p0]];
3927 MInt vB2 = edges[startSetIndex][prevEdge].vertices[id_p0];
3928 MInt vA = edges[startSetIndex][nextEdge].vertices[id_p1];
3929 MInt vB = edges[startSetIndex][nextEdge].vertices[otherId[id_p1]];
3932 MInt vO = otherVert;
3941 MFloat vec1[nDim] = {F0, F0, F0};
3942 MFloat vec2[nDim] = {F0, F0, F0};
3943 MFloat vec3[nDim] = {F0, F0, F0};
3944 MFloat ori[nDim] = {F0, F0, F0};
3945 MFloat normal[3] = {F0, F0, F0};
3946 MFloat area0 = faces[startSetIndex][face0].area;
3947 MFloat area1 = faces[startSetIndex][face1].area;
3948 ASSERT(area0 + area1 > F0,
"");
3949 for(
MInt i = 0; i < nDim; i++) {
3950 normal[i] = (area0 * faces[startSetIndex][face0].normal[i]
3951 + area1 * faces[startSetIndex][face1].normal[i])
3954 for(
MInt i = 0; i < nDim; i++) {
3956 vertices[startSetIndex][vB2].coordinates[i] - vertices[startSetIndex][vA2].coordinates[i];
3958 vertices[startSetIndex][vB].coordinates[i] - vertices[startSetIndex][vA].coordinates[i];
3960 vertices[startSetIndex][vO].coordinates[i] - vertices[startSetIndex][vV].coordinates[i];
3972 dotp /= (abs1 * abs2);
3973 for(
MInt i = 0; i < nDim; i++) {
3978 crossProduct(ori, vec1, vec2);
3979 for(
MInt i = 0; i < nDim; i++) {
3980 dotp2 += ori[i] * normal[i];
3982 crossProduct(ori, vec1, vec3);
3983 for(
MInt i = 0; i < nDim; i++) {
3984 dotp3 += ori[i] * normal[i];
3986 crossProduct(ori, vec2, vec3);
3987 for(
MInt i = 0; i < nDim; i++) {
3988 dotp4 += ori[i] * normal[i];
3991 if(dotp2 < -0.1 && fabs(F1 - dotp) > 1e-2
3992 &&
mMin(area0, area1) > 1e-5 * pow(cellLength0, (
MFloat)(nDim - 1))) {
3993 if(noLocBodySrfcs + concaveCnt > maxNoSurfaces) {
3994 cerr <<
"Warning: removal of concave polygons was skipped, since "
3995 "FvBndryCell<nDim>::m_maxNoSurfaces is not large enough."
4000 if(dotp3 > -1e-3 && dotp4 > -1e-3) {
4001 auto it2 = pureBodyEdges.begin();
4002 while(it2 != pureBodyEdges.end()) {
4005 MInt bedge = (*it2);
4006 if(edgeCutCellPointer[bedge] != c) {
4010 MInt faceA = edges[startSetIndex][bedge].face[0];
4011 MInt faceB = edges[startSetIndex][bedge].face[1];
4012 if((face0 == faceA && face1 == faceB) || (face1 == faceA && face0 == faceB)) {
4013 it2 = pureBodyEdges.erase(it2);
4031 MBool somethingChanged =
true;
4032 MInt currentMultiFace = 0;
4033 while(!pureBodyEdges.empty()) {
4034 MInt pbedge = pureBodyEdges.front();
4035 pureBodyEdges.pop_front();
4036 multiFaces.emplace_back();
4037 MInt pbface0 = edges[startSetIndex][pbedge].face[0];
4038 multiFaceConnection[pbface0] = currentMultiFace;
4039 multiFaces[currentMultiFace].faces.push_back(pbface0);
4040 multiFaces[currentMultiFace].bodyId = faces[startSetIndex][pbface0].bodyId;
4041 MInt pbface1 = edges[startSetIndex][pbedge].face[1];
4042 multiFaceConnection[pbface1] = currentMultiFace;
4043 multiFaces[currentMultiFace].faces.push_back(pbface1);
4044 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][pbface0].edges.size(); e++) {
4045 if(faces[startSetIndex][pbface0].edges[e].first != pbedge) {
4046 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface0].edges[e].first,
4047 faces[startSetIndex][pbface0].edges[e].second);
4051 for(
MInt ee = 0; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4052 if(faces[startSetIndex][pbface1].edges[ee].first == pbedge) {
4057 for(
MInt ee = 1; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4058 MInt eee = (edgeIndex + ee) % faces[startSetIndex][pbface1].edges.
size();
4059 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface1].edges[eee].first,
4060 faces[startSetIndex][pbface1].edges[eee].second);
4064 somethingChanged =
true;
4066 while(somethingChanged) {
4067 somethingChanged =
false;
4068 auto it = pureBodyEdges.begin();
4069 while(it != pureBodyEdges.end()) {
4072 if(edgeCutCellPointer[edge] != c) {
4076 MInt face0 = edges[startSetIndex][edge].face[0];
4077 MInt face1 = edges[startSetIndex][edge].face[1];
4078 if(multiFaceConnection[face0] == currentMultiFace && multiFaceConnection[face1] == currentMultiFace) {
4079 it = pureBodyEdges.erase(it);
4081 if(multiFaceConnection[face0] == currentMultiFace) {
4082 multiFaceConnection[face1] = currentMultiFace;
4083 it = pureBodyEdges.erase(it);
4084 multiFaces[currentMultiFace].faces.push_back(face1);
4085 auto it2 = multiFaces[currentMultiFace].edges.
begin();
4086 for(it2 = multiFaces[currentMultiFace].edges.begin();
4087 it2 != multiFaces[currentMultiFace].edges.end();
4089 MInt edge2 = (*it2).first;
4096 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4097 if(faces[startSetIndex][face1].edges[e].first == edge) {
4102 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4103 MInt ee = (edgeIndex + e) % faces[startSetIndex][face1].edges.size();
4104 multiFaces[currentMultiFace].edges.insert(
4105 it2, make_pair(faces[startSetIndex][face1].edges[ee].first,
4106 faces[startSetIndex][face1].edges[ee].second));
4108 multiFaces[currentMultiFace].edges.erase(it2);
4109 somethingChanged =
true;
4110 }
else if(multiFaceConnection[face1] == currentMultiFace) {
4111 multiFaceConnection[face0] = currentMultiFace;
4112 it = pureBodyEdges.erase(it);
4113 multiFaces[currentMultiFace].faces.push_back(face0);
4114 auto it2 = multiFaces[currentMultiFace].edges.
begin();
4115 for(it2 = multiFaces[currentMultiFace].edges.begin();
4116 it2 != multiFaces[currentMultiFace].edges.end();
4118 MInt edge2 = (*it2).first;
4125 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4126 if(faces[startSetIndex][face0].edges[e].first == edge) {
4131 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4132 MInt ee = (edgeIndex + e) % faces[startSetIndex][face0].edges.size();
4133 multiFaces[currentMultiFace].edges.insert(
4134 it2, make_pair(faces[startSetIndex][face0].edges[ee].first,
4135 faces[startSetIndex][face0].edges[ee].second));
4137 multiFaces[currentMultiFace].edges.erase(it2);
4138 somethingChanged =
true;
4150 for(
MInt e = 0; (unsigned)e < edges[startSetIndex].size(); e++) {
4151 if(edgeCutCellPointer[e] != c)
continue;
4152 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 1
4153 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 1) {
4154 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].bodyId
4155 == faces[startSetIndex][edges[startSetIndex][e].face[1]].bodyId) {
4156 pureBodyEdges.push_back(e);
4158 }
else if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceType == 0
4159 && faces[startSetIndex][edges[startSetIndex][e].face[1]].faceType == 0) {
4160 if(faces[startSetIndex][edges[startSetIndex][e].face[0]].faceId
4161 == faces[startSetIndex][edges[startSetIndex][e].face[1]].faceId) {
4162 pureBodyEdges.push_back(e);
4167 const MInt faceNum = faces[startSetIndex].size();
4168 ASSERT(faceNum <= maxNoFaces,
"");
4169 for(
MInt i = 0; i < faceNum; i++)
4170 for(
MInt j = 0; j < faceNum; j++)
4171 normalDotProduct(i, j) = F0;
4172 MInt noBodyFaces = 0;
4173 for(
MInt f = 0; (unsigned)f < cutCells[c].faces.size(); f++) {
4174 MInt face = cutCells[c].faces[f];
4175 if(faces[startSetIndex][face].faceType == 1) bodyFaces[noBodyFaces++] = face;
4177 for(
MInt f = 0; f < noBodyFaces; f++) {
4178 MInt face1 = bodyFaces[f];
4179 for(
MInt ff = 0; ff <= f; ff++) {
4180 MInt face2 = bodyFaces[ff];
4181 for(
MInt i = 0; i < nDim; i++)
4182 normalDotProduct(face1, face2) +=
4183 faces[startSetIndex][face1].normal[i] * faces[startSetIndex][face2].normal[i];
4184 normalDotProduct(face1, face2) = F1 - normalDotProduct(face1, face2);
4185 normalDotProduct(face2, face1) = normalDotProduct(face1, face2);
4192 MBool somethingChanged =
true;
4193 MInt currentMultiFace = 0;
4194 while(!pureBodyEdges.empty()) {
4195 MInt pbedge = pureBodyEdges.front();
4196 pureBodyEdges.pop_front();
4197 MInt pbface0 = edges[startSetIndex][pbedge].face[0];
4198 MInt pbface1 = edges[startSetIndex][pbedge].face[1];
4199 if(normalDotProduct(pbface0, pbface1) > m_bodyFaceJoinCriterion)
continue;
4200 multiFaces.emplace_back();
4201 multiFaceConnection[pbface0] = currentMultiFace;
4202 multiFaces[currentMultiFace].faces.push_back(pbface0);
4203 multiFaces[currentMultiFace].bodyId = faces[startSetIndex][pbface0].bodyId;
4204 multiFaceConnection[pbface1] = currentMultiFace;
4205 multiFaces[currentMultiFace].faces.push_back(pbface1);
4206 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][pbface0].edges.size(); e++) {
4207 if(faces[startSetIndex][pbface0].edges[e].first != pbedge) {
4208 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface0].edges[e].first,
4209 faces[startSetIndex][pbface0].edges[e].second);
4213 for(
MInt ee = 0; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4214 if(faces[startSetIndex][pbface1].edges[ee].first == pbedge) {
4219 for(
MInt ee = 1; (unsigned)ee < faces[startSetIndex][pbface1].edges.size(); ee++) {
4220 MInt eee = (edgeIndex + ee) % faces[startSetIndex][pbface1].edges.
size();
4221 multiFaces[currentMultiFace].edges.emplace_back(faces[startSetIndex][pbface1].edges[eee].first,
4222 faces[startSetIndex][pbface1].edges[eee].second);
4226 somethingChanged =
true;
4228 while(somethingChanged) {
4229 somethingChanged =
false;
4230 auto it = pureBodyEdges.begin();
4231 while(it != pureBodyEdges.end()) {
4234 if(edgeCutCellPointer[edge] != c) {
4238 MInt face0 = edges[startSetIndex][edge].face[0];
4239 MInt face1 = edges[startSetIndex][edge].face[1];
4240 if(multiFaceConnection[face0] == currentMultiFace && multiFaceConnection[face1] == currentMultiFace) {
4241 it = pureBodyEdges.erase(it);
4242 }
else if(multiFaceConnection[face0] == currentMultiFace) {
4243 it = pureBodyEdges.erase(it);
4244 MBool mayBeConnected =
true;
4245 for(
MInt nf = 0; (unsigned)nf < multiFaces[currentMultiFace].faces.size(); nf++) {
4246 MInt nFace = multiFaces[currentMultiFace].faces[nf];
4247 if(normalDotProduct(face1, nFace) > m_bodyFaceJoinCriterion) {
4248 mayBeConnected =
false;
4252 if(!mayBeConnected)
continue;
4253 multiFaceConnection[face1] = currentMultiFace;
4254 multiFaces[currentMultiFace].faces.push_back(face1);
4255 auto it2 = multiFaces[currentMultiFace].edges.
begin();
4256 for(it2 = multiFaces[currentMultiFace].edges.begin();
4257 it2 != multiFaces[currentMultiFace].edges.end();
4259 MInt edge2 = (*it2).first;
4266 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4267 if(faces[startSetIndex][face1].edges[e].first == edge) {
4272 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face1].edges.size(); e++) {
4273 MInt ee = (edgeIndex + e) % faces[startSetIndex][face1].edges.size();
4274 multiFaces[currentMultiFace].edges.insert(
4275 it2, make_pair(faces[startSetIndex][face1].edges[ee].first,
4276 faces[startSetIndex][face1].edges[ee].second));
4278 multiFaces[currentMultiFace].edges.erase(it2);
4279 somethingChanged =
true;
4280 }
else if(multiFaceConnection[face1] == currentMultiFace) {
4281 it = pureBodyEdges.erase(it);
4282 MBool mayBeConnected =
true;
4283 for(
MInt nf = 0; (unsigned)nf < multiFaces[currentMultiFace].faces.size(); nf++) {
4284 MInt nFace = multiFaces[currentMultiFace].faces[nf];
4285 if(normalDotProduct(face0, nFace) > m_bodyFaceJoinCriterion) {
4286 mayBeConnected =
false;
4290 if(!mayBeConnected)
continue;
4291 multiFaceConnection[face0] = currentMultiFace;
4292 multiFaces[currentMultiFace].faces.push_back(face0);
4293 auto it2 = multiFaces[currentMultiFace].edges.
begin();
4294 for(it2 = multiFaces[currentMultiFace].edges.begin();
4295 it2 != multiFaces[currentMultiFace].edges.end();
4297 MInt edge2 = (*it2).first;
4304 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4305 if(faces[startSetIndex][face0].edges[e].first == edge) {
4310 for(
MInt e = 1; (unsigned)e < faces[startSetIndex][face0].edges.size(); e++) {
4311 MInt ee = (edgeIndex + e) % faces[startSetIndex][face0].edges.size();
4312 multiFaces[currentMultiFace].edges.insert(
4313 it2, make_pair(faces[startSetIndex][face0].edges[ee].first,
4314 faces[startSetIndex][face0].edges[ee].second));
4316 multiFaces[currentMultiFace].edges.erase(it2);
4317 somethingChanged =
true;
4328 mTerm(1, AT_,
"ERROR: invalid bodyFaceJoinMode specified. exiting...");
4333 for(
MInt mf = 0; (unsigned)mf < multiFaces.size(); mf++) {
4334 MFloat normal[3] = {F0, F0, F0};
4335 MFloat center[3] = {F0, F0, F0};
4337 for(
MInt f = 0; (unsigned)f < multiFaces[mf].faces.size(); f++) {
4338 const MInt face = multiFaces[mf].faces[f];
4339 MFloat faceArea = faces[startSetIndex][face].area;
4341 if(faceArea < m_eps) {
4345 for(
MInt i = 0; i < nDim; i++) {
4346 normal[i] += faces[startSetIndex][face].normal[i] * faceArea;
4347 center[i] += faces[startSetIndex][face].center[i] * faceArea;
4351 for(
MInt i = 0; i < nDim; i++) {
4352 absNorm += normal[i] * normal[i];
4354 absNorm = sqrt(absNorm);
4355 for(
MInt i = 0; i < nDim; i++) {
4356 multiFaces[mf].center[i] = center[i] / area2;
4357 ASSERT(!std::isnan(multiFaces[mf].center[i]),
"");
4358 multiFaces[mf].normal[i] = normal[i] / absNorm;
4359 ASSERT(!std::isnan(multiFaces[mf].normal[i]),
"");
4361 multiFaces[mf].area = absNorm;
4365 if(!multiFaces.empty()) {
4368 for(
MInt f = 0; (unsigned)f < cutCells[c].faces.size(); f++) {
4369 MInt face = cutCells[c].faces[f];
4370 if(multiFaceConnection[face] == -1) tmp_faces[noFaces++] = cutCells[c].faces[f];
4372 cutCells[c].faces.clear();
4373 for(
MInt f = 0; f < noFaces; f++) {
4374 cutCells[c].faces.push_back(tmp_faces[f]);
4379 for(
MInt mf = 0; (unsigned)mf < multiFaces.size(); mf++) {
4382 const MInt newFace = faces[startSetIndex].
size();
4383 const MInt faceId = faces[startSetIndex][multiFaces[mf].faces[0]].faceId;
4384 const MInt faceType = faces[startSetIndex][multiFaces[mf].faces[0]].faceType;
4385 const MInt bodyId = faces[startSetIndex][multiFaces[mf].faces[0]].bodyId;
4387 faces[startSetIndex].emplace_back(faceId, faceType, bodyId);
4388 faces[startSetIndex][newFace].area = multiFaces[mf].area;
4389 for(
MInt i = 0; i < nDim; i++) {
4390 faces[startSetIndex][newFace].center[i] = multiFaces[mf].center[i];
4391 faces[startSetIndex][newFace].normal[i] = multiFaces[mf].normal[i];
4393 faces[startSetIndex][newFace].w = F0;
4394 faces[startSetIndex][newFace].cutCell = c;
4398 MBool somethingChanged =
true;
4399 while(somethingChanged) {
4400 somethingChanged =
false;
4403 auto it = multiFaces[mf].edges.begin();
4404 while(it != multiFaces[mf].edges.end()) {
4405 auto itLast = multiFaces[mf].edges.end();
4406 if(it == multiFaces[mf].edges.begin()) {
4407 itLast = multiFaces[mf].edges.end();
4412 ASSERT(itLast != it,
"");
4413 if((it->first == itLast->first) && (it->second == -(itLast->second))) {
4414 multiFaces[mf].edges.erase(itLast);
4415 it = multiFaces[mf].edges.erase(it);
4416 somethingChanged =
true;
4425 for(
auto& edge : multiFaces[mf].edges) {
4426 faces[startSetIndex][newFace].edges.emplace_back(edge.first, edge.second);
4427 MInt vertexId = edges[startSetIndex][edge.first].vertices[0];
4428 if(edge.second == -1) vertexId = edges[startSetIndex][edge.first].vertices[1];
4429 faces[startSetIndex][newFace].vertices.push_back(vertexId);
4431 cutCells[c].faces.push_back(newFace);
4435 RECORD_TIMER_STOP(tCutFace_5b);
4444 for(
MInt sc = 0; sc < (
MInt)cutCells.size(); sc++) {
4445 if(!cutCells[sc].faces.empty()) {
4448 MInt vertexRemap[maxNoVertices];
4449 fill_n(vertexRemap, maxNoVertices, -1);
4450 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4451 MInt face = cutCells[sc].faces[f];
4452 noEdges += faces[startSetIndex][face].vertices.size();
4453 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
4454 const MInt vertex = faces[startSetIndex][face].vertices[e];
4455 if(vertexRemap[vertex] < 0) {
4456 vertexRemap[vertex] = pcnt;
4462 if(pcnt + cutCells[sc].faces.size() != 2 + noEdges / 2 || noEdges % 2 == 1) {
4463 cerr << grid().domainId() <<
"Euler's polyhedral formula (1) not fulfilled ("
4464 << pcnt + (signed)cutCells[sc].faces.size() - 2 - noEdges / 2 <<
") " << pcnt <<
" "
4465 << cutCells[sc].faces.size() <<
" " << noEdges / 2 <<
" " << noEdges <<
" for cell " << cellId <<
" "
4466 <<
globalTimeStep <<
" " << faces[startSetIndex].size() <<
" " << cutCells.size() <<
" "
4467 << grid().tree().globalId(cellId) <<
" " << sc << endl;
4468 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], -1);
4471 cerr <<
"Deleting splitchilds " << endl;
4472 cutCells.erase(cutCells.begin() + sc);
4487 MInt splitCellId = -1;
4488 MInt noSplitChildren = cutCells.size();
4489 if(noSplitChildren > CC::maxSplitCells) {
4490 mTerm(1, AT_,
"Too many split cells generated, see CutCell::maxSplitCells");
4492 cutCellData[cutc].noSplitChilds = 0;
4495 if(noSplitChildren > 1) {
4500 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
4501 MUint splitcutc = cutCellData.size();
4502 cutCellData[cutc].splitChildIds[cutCellData[cutc].noSplitChilds++] = splitcutc;
4503 cutCellData.emplace_back();
4506 cutCellData[splitcutc].cellId = -1;
4509 cutCellData[splitcutc].splitParentId = cutc;
4540 ASSERT(cutCellData[cutc].noSplitChilds == 0 || cutCellData[cutc].noSplitChilds > 1,
"");
4542 if(noSplitChildren > 1) {
4543 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
4544 splitCellList[sc] = cutCellData[cutc].splitChildIds[sc];
4547 splitCellList[0] = cutc;
4552 MInt noFacesPerCartesianFaceAll[6] = {0, 0, 0, 0, 0, 0};
4553 MInt noBodySurfacesAll = 0;
4557 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
4558 RECORD_TIMER_START(tCutFace_6);
4570 MInt cutCellId = splitCellList[sc];
4576 if(cutCellData[cutCellId].splitParentId > -1) {
4577 for(
MInt corner = 0; corner < m_noCorners; corner++) {
4578 cutCellData[cutCellId].cornerIsInsideGeometry[0][corner] =
true;
4580 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4581 MInt face = cutCells[sc].faces[f];
4582 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
4583 const MInt vStart = faces[startSetIndex][face].vertices[e];
4584 if(vertices[startSetIndex][vStart].
pointType == 0) {
4585 MInt corner = vertices[startSetIndex][vStart].pointId;
4586 cutCellData[cutCellId].cornerIsInsideGeometry[0][corner] =
false;
4595 MInt noFacesPerCartesianFace[6] = {0, 0, 0, 0, 0, 0};
4596 MInt noBodySurfaces = 0;
4597 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4598 MInt face = cutCells[sc].faces[f];
4599 if(faces[startSetIndex][face].faceType == 0) {
4600 noFacesPerCartesianFace[faces[startSetIndex][face].faceId]++;
4601 noFacesPerCartesianFaceAll[faces[startSetIndex][face].faceId]++;
4604 noBodySurfacesAll++;
4647 if(noBodySurfaces > maxNoSurfaces || noBodySurfaces > CC::maxNoBoundarySurfaces) {
4648 cerr <<
"more than FvBndryCell<nDim>::m_maxNoSurfaces cut surfaces detected for a cell. This is not "
4649 "implemented in MB framework yet. cell "
4650 << cellId <<
" " << grid().domainId() <<
" " <<
globalTimeStep <<
", splitChild " << sc <<
" "
4651 << noBodySurfaces <<
" maximum " << maxNoSurfaces <<
" " << CC::maxNoBoundarySurfaces <<
" "
4652 << noSplitChildren <<
" " << faces[startSetIndex].size() << endl;
4653 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], startSetIndex);
4655 " more than 3 cut surfaces detected for a cell. This is not implemented in MB framework yet. "
4664 if(cutCells[sc].faces.empty()) {
4672 cutCellData[cutCellId].volume = F0;
4676 if(cutCells[sc].volume > a_gridCellVolume(cutCellId)) {
4677 cutCells[sc].volume = a_gridCellVolume(cutCellId);
4683 for(
MInt f = 0; f < m_noDirs; f++) {
4684 if(noFacesPerCartesianFace[f] == 0) {
4686 cutCellData[cutCellId].externalFaces[f] =
true;
4690 if(splitCellId < 0) {
4730 ASSERT(vertices[startSetIndex].size() <= (
unsigned)maxNoVertices,
"");
4731 for(
MInt cp = 0; cp < maxNoVertices; cp++)
4732 newCutPointId[cp] = -1;
4733 for(
MInt i = 0; (unsigned)i < vertices[startSetIndex].size(); i++)
4734 isPartOfBodySurface[i] =
false;
4735 for(
MInt e = 0; e < 2 * m_noEdges; e++) {
4737 pointEdgeId[e] = cutCellData[cutc].cutEdges[e];
4741 ASSERT(cutCells[sc].volume >= F0,
"");
4742 ASSERT(!(std::isnan(cutCells[sc].volume)),
"");
4743 ASSERT(!(std::isnan(cutCells[sc].center[sc])),
"");
4744 ASSERT(!(std::isnan(cutCells[sc].center[1])),
"");
4745 ASSERT(!(std::isnan(cutCells[sc].center[2])),
"");
4746 ASSERT(!(std::isinf(cutCells[sc].volume)),
"");
4747 ASSERT(!(std::isinf(cutCells[sc].center[sc])),
"");
4748 ASSERT(!(std::isinf(cutCells[sc].center[1])),
"");
4749 ASSERT(!(std::isinf(cutCells[sc].center[2])),
"");
4752 cutCellData[cutCellId].volume = cutCells[sc].volume;
4758 for(
MInt i = 0; i < nDim; i++) {
4760 cutCellData[cutCellId].volumetricCentroid[i] = cutCells[sc].center[i] - a_coordinate(cellId, i);
4764 cutCellData[cutCellId].noCartesianSurfaces = 0;
4765 cutCellData[cutCellId].noBoundarySurfaces = 0;
4769 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
4770 MInt face = cutCells[sc].faces[f];
4771 ASSERT(faces[startSetIndex][face].area >= F0,
"");
4772 ASSERT(!(std::isnan(faces[startSetIndex][face].area)),
"");
4773 ASSERT(!(std::isnan(faces[startSetIndex][face].center[0])),
"");
4774 ASSERT(!(std::isnan(faces[startSetIndex][face].center[1])),
"");
4775 ASSERT(!(std::isnan(faces[startSetIndex][face].center[2])),
"");
4776 ASSERT(!(std::isinf(faces[startSetIndex][face].area)),
"");
4777 ASSERT(!(std::isinf(faces[startSetIndex][face].center[0])),
"");
4778 ASSERT(!(std::isinf(faces[startSetIndex][face].center[1])),
"");
4779 ASSERT(!(std::isinf(faces[startSetIndex][face].center[2])),
"");
4785 if(faces[startSetIndex][face].faceType == 1) {
4788 cutCellData[cutCellId].boundarySurfaceArea[noBodySurfaces] = faces[startSetIndex][face].area;
4790 for(
MInt i = 0; i < nDim; i++) {
4793 cutCellData[cutCellId].boundarySurfaceCentroid[noBodySurfaces][i] = faces[startSetIndex][face].center[i];
4794 cutCellData[cutCellId].boundarySurfaceNormal[noBodySurfaces][i] = faces[startSetIndex][face].normal[i];
4799 cutCellData[cutCellId].boundarySurfaceBndryCndId[noBodySurfaces] = -1;
4800 cutCellData[cutCellId].boundarySurfaceBodyId[noBodySurfaces] = faces[startSetIndex][face].bodyId;
4803 for(
MInt v = 0; (unsigned)v < faces[startSetIndex][face].vertices.size(); v++) {
4804 const MInt vertex = faces[startSetIndex][face].vertices[v];
4806 isPartOfBodySurface[vertex] =
true;
4808 newCutPointId[vertex] = v;
4809 vertices[startSetIndex][vertex].cartSrfcId = noBodySurfaces;
4836 cutCellData[cutCellId].noBoundarySurfaces++;
4837 ASSERT(cutCellData[cutCellId].noBoundarySurfaces == noBodySurfaces,
"");
4844 const MInt dirId = faces[startSetIndex][face].faceId;
4845 const MInt spaceId = dirId / 2;
4902 if(dirId % 2 == 0) sign = -F1;
4905 for(
MInt i = 0; i < nDim; i++) {
4908 cutCellData[cutCellId].cartFaceCentroid[cutCellData[cutCellId].noCartesianSurfaces][i] =
4909 faces[startSetIndex][face].center[i];
4913 cutCellData[cutCellId].cartFaceCentroid[cutCellData[cutCellId].noCartesianSurfaces][spaceId] =
4914 a_coordinate(cellId, spaceId) + sign * cellHalfLength;
4917 cutCellData[cutCellId].cartFaceArea[cutCellData[cutCellId].noCartesianSurfaces] =
4918 faces[startSetIndex][face].area;
4919 cutCellData[cutCellId].cartFaceDir[cutCellData[cutCellId].noCartesianSurfaces] =
4920 faces[startSetIndex][face].faceId;
4981 cutCellData[cutCellId].noCartesianSurfaces++;
4982 if(cutCellData[cutCellId].noCartesianSurfaces > CC::maxNoCartesianSurfaces) {
4983 mTerm(1, AT_,
"Increase CutCell::maxNoCartesianSurfaces");
4988 if(noSplitChildren > 1) {
4990 cutCellData[cutc].noBoundarySurfaces = 0;
4993 cutCellData[cutCellId].noBoundarySurfaces = noBodySurfaces;
4995 RECORD_TIMER_STOP(tCutFace_6);
4997 for(
MInt k = 0; k < CC::noFaces; k++) {
4998 cutCellData[cutCellId].noFacesPerCartesianDir[k] = noFacesPerCartesianFace[k];
5003 RECORD_TIMER_START(tCutFace_7);
5007 cutCellData[cutCellId].noTotalFaces = 0;
5008 cutCellData[cutCellId].noAdditionalVertices = 0;
5009 MInt vertexMap[maxNoVertices];
5010 fill_n(vertexMap, maxNoVertices, -1);
5012 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
5013 MInt face = cutCells[sc].faces[f];
5016 const MInt facecnt = cutCellData[cutCellId].noTotalFaces;
5017 cutCellData[cutCellId].allFacesNoPoints[facecnt] = 0;
5018 cutCellData[cutCellId].allFacesBodyId[facecnt] = faces[startSetIndex][face].bodyId;
5020 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
5021 const MInt vertex = faces[startSetIndex][face].vertices[e];
5022 MInt pointId = vertices[startSetIndex][vertex].pointId;
5023 MInt storePointId = -1;
5024 if(vertices[startSetIndex][vertex].
pointType == 0) {
5028 ASSERT(pointId < CC::noCorners,
"");
5029 storePointId = pointId;
5030 }
else if(pointId > -1) {
5031 ASSERT(vertices[startSetIndex][vertex].
pointType == 1,
"");
5043 ASSERT(pointId < cutCellData[cutc].noCutPoints,
"");
5044 storePointId = CC::noCorners + pointId;
5046 ASSERT(vertices[startSetIndex][vertex].
pointType > 1,
"");
5047 if(vertexMap[vertex] < 0) {
5048 for(
MInt i = 0; i < nDim; i++) {
5049 cutCellData[cutCellId].additionalVertices[cutCellData[cutCellId].noAdditionalVertices][i] =
5050 vertices[startSetIndex][vertex].coordinates[i];
5052 vertexMap[vertex] = cutCellData[cutCellId].noAdditionalVertices;
5053 cutCellData[cutCellId].noAdditionalVertices++;
5054 if(cutCellData[cutCellId].noAdditionalVertices > CC::maxNoAdditionalVertices) {
5055 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex],
globalTimeStep);
5057 ASSERT(cutCellData[cutCellId].noAdditionalVertices <= CC::maxNoAdditionalVertices,
"");
5059 storePointId = CC::noCorners + CC::maxNoCutPoints + vertexMap[vertex];
5061 cutCellData[cutCellId].allFacesPointIds[facecnt][cutCellData[cutCellId].allFacesNoPoints[facecnt]] =
5063 cutCellData[cutCellId].allFacesNoPoints[facecnt]++;
5064 ASSERT(cutCellData[cutCellId].allFacesNoPoints[facecnt] <= CC::maxNoFaceVertices,
5065 "has: " + to_string(cutCellData[cutCellId].allFacesNoPoints[facecnt]) +
" max is "
5066 + to_string(CC::maxNoFaceVertices));
5069 cutCellData[cutCellId].noTotalFaces++;
5198 if(!cutCells[sc].faces.empty()) {
5201 MInt vertexRemap[maxNoVertices];
5202 fill_n(vertexRemap, maxNoVertices, -1);
5203 for(
MInt f = 0; (unsigned)f < cutCells[sc].faces.size(); f++) {
5204 MInt face = cutCells[sc].faces[f];
5206 noEdges += faces[startSetIndex][face].vertices.size();
5207 for(
MInt e = 0; (unsigned)e < faces[startSetIndex][face].vertices.size(); e++) {
5208 const MInt vertex = faces[startSetIndex][face].vertices[e];
5209 if(vertexRemap[vertex] < 0) {
5210 vertexRemap[vertex] = pcnt;
5215 if(pcnt + cutCells[sc].faces.size() != 2 + noEdges / 2 || noEdges % 2 == 1) {
5216 cerr << grid().domainId() <<
"Euler's polyhedral formula not fulfilled ("
5217 << pcnt + (signed)cutCells[sc].faces.size() - 2 - noEdges / 2 <<
") " << pcnt <<
" "
5218 << cutCells[sc].faces.size() <<
" " << noEdges / 2 <<
" " << noEdges <<
" for cell " << cellId <<
" "
5219 <<
globalTimeStep <<
" " << faces[startSetIndex].size() <<
" " << cutCells.size() <<
" "
5220 << grid().tree().globalId(cellId) <<
" " << sc << endl;
5221 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], -1);
5224 RECORD_TIMER_STOP(tCutFace_7);
5228 if(noSplitChildren > 1) {
5229 cutCellData[cutc].volume = F0;
5230 for(
MInt sc = 0; sc < noSplitChildren; sc++) {
5231 MInt cutCellId = splitCellList[sc];
5232 cutCellData[cutc].volume += cutCellData[cutCellId].volume;
5234 for(
MInt dir = 0; dir < m_noDirs; dir++) {
5235 cutCellData[cutc].noFacesPerCartesianDir[dir] = noFacesPerCartesianFaceAll[dir];
5237 if(cutCellData[cutc].noFacesPerCartesianDir[dir] == 0) {
5238 cutCellData[cutc].externalFaces[dir] =
true;
5453 if(
globalTimeStep == debugTimeStep && cellId == debugCellId && grid().domainId() == debugDomainId) {
5454 writeVTKFileOfCell(cellId, &faces[startSetIndex], &vertices[startSetIndex], -2);
5456 const MChar* fileName =
"cell-Info_";
5457 stringstream fileName2;
5458 fileName2 << fileName << grid().tree().globalId(cellId) <<
".txt";
5461 ofl.open((fileName2.str()).c_str(), ofstream::trunc);
5464 ofl.setf(ios::fixed);
5467 ofl <<
"Writing Cell " << grid().tree().globalId(cellId) <<
" " << cellId <<
" " << cutc <<
" noCuts "
5468 << noIndividualCuts <<
" noCutCells " << cutCells.size() <<
" vol " << cutCellData[cutc].volume
5469 <<
" SplitChilds " << cutCellData[cutc].noSplitChilds <<
" " << cutCellData[cutc].splitParentId << endl
5471 for(
MInt dir = 0; dir < m_noDirs; dir++) {
5472 ofl <<
"Dir " << dir <<
" external " << cutCellData[cutc].externalFaces[dir] <<
" faces "
5473 << cutCellData[cutc].noFacesPerCartesianDir[dir] << endl
5477 ofl <<
" Cartesian-Surf: " << cutCellData[cutc].noCartesianSurfaces << endl << endl;
5479 for(
MInt id = 0;
id < cutCellData[cutc].noCartesianSurfaces;
id++) {
5480 ofl <<
"Id " <<
id <<
" dir " << cutCellData[cutc].cartFaceDir[
id] <<
" area "
5481 << cutCellData[cutc].cartFaceArea[
id] <<
" coord " << cutCellData[cutc].cartFaceCentroid[
id][0]
5482 << cutCellData[cutc].cartFaceCentroid[
id][1] << cutCellData[cutc].cartFaceCentroid[
id][2] << endl
5486 ofl <<
" Bndry-Surfaces: " << cutCellData[cutc].noBoundarySurfaces << endl << endl;
5488 for(
MInt id = 0;
id < cutCellData[cutc].noBoundarySurfaces;
id++) {
5489 ofl <<
"Id " <<
id <<
" normal " << cutCellData[cutc].boundarySurfaceNormal[
id][0]
5490 << cutCellData[cutc].boundarySurfaceNormal[
id][1] << cutCellData[cutc].boundarySurfaceNormal[
id][2]
5491 <<
" center " << cutCellData[cutc].boundarySurfaceCentroid[
id][0]
5492 << cutCellData[cutc].boundarySurfaceCentroid[
id][1] << cutCellData[cutc].boundarySurfaceCentroid[
id][2]
5493 <<
" area " << cutCellData[cutc].boundarySurfaceArea[
id] << endl
5533template <MInt nDim_>
5535 const MInt* candidateIds,
5536 std::vector<MInt>& candidatesOrder) {
5539 const MInt edgesPerDim =
IPOW2(nDim_ - 1);
5541 const MInt DOFStencil[12] = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
5544 const MInt faceStencil[2][12] = {{0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1}, {4, 4, 4, 4, 5, 5, 5, 5, 2, 2, 3, 3}};
5551 const MInt reverseEdgeStencil[3][12] = {{1, 0, 3, 2, 5, 4, 7, 6, 9, 8, 11, 10},
5552 {4, 5, 6, 7, 0, 1, 2, 3, 10, 11, 8, 9},
5553 {5, 4, 7, 6, 1, 0, 3, 2, 11, 10, 9, 8}};
5555 const MInt signStencil[3][12] = {{-1, 1, 0, 0, -1, 1, 0, 0, -1, 1, -1, 1},
5556 {0, 0, -1, 1, 0, 0, -1, 1, -1, -1, 1, 1},
5557 {-1, -1, -1, -1, 1, 1, 1, 1, 0, 0, 0, 0}};
5561 const MInt nodeStencil[2][12] = {{0, 1, 0, 2, 4, 5, 4, 6, 0, 1, 2, 3}, {2, 3, 1, 3, 6, 7, 5, 7, 4, 5, 6, 7}};
5566 const MFloat eps = m_multiCutCell ? 10 * m_eps : m_eps;
5573 for(
MInt cndt = 0; cndt < (signed)candidates.size(); cndt++) {
5574 const MInt cellId = candidates[cndt].cellId;
5575 const MFloat cellLength = a_cellLengthAtCell(cellId);
5576 const MFloat cellHalfLength = F1B2 * cellLength;
5577 for(
MInt edge = 0; edge < m_noEdges; edge++) {
5578 if(candidates[cndt].edgeChecked[edge])
continue;
5581 MInt face[2]{-1, -1};
5583 face[0] = faceStencil[0][edge];
5585 IF_CONSTEXPR(nDim_ == 3) {
5587 face[1] = faceStencil[1][edge];
5591 const MInt edgeDOF = DOFStencil[edge];
5594 MInt directNeighbourIds[4]{-1, -1, -1, -1};
5595 directNeighbourIds[0] = cellId;
5597 if(grid().tree().hasNeighbor(cellId, face[0]) > 0) {
5598 directNeighbourIds[1] = grid().tree().neighbor(cellId, face[0]);
5600 MInt reverseEdge[edgesPerDim];
5601 reverseEdge[0] = edge;
5602 reverseEdge[1] = reverseEdgeStencil[0][edge];
5604 IF_CONSTEXPR(nDim_ == 3) {
5605 if(grid().tree().hasNeighbor(cellId, face[1]) > 0) {
5606 directNeighbourIds[2] = grid().tree().neighbor(cellId, face[1]);
5609 if(grid().tree().hasNeighbor(cellId, face[0]) > 0) {
5610 if(grid().tree().hasNeighbor(grid().tree().neighbor(cellId, face[0]), face[1]) > 0) {
5611 directNeighbourIds[3] = grid().tree().neighbor(grid().tree().neighbor(cellId, face[0]), face[1]);
5614 if(directNeighbourIds[3] < 0) {
5615 if(grid().tree().hasNeighbor(cellId, face[1]) > 0) {
5616 if(grid().tree().hasNeighbor(grid().tree().neighbor(cellId, face[1]), face[0]) > 0) {
5617 directNeighbourIds[3] = grid().tree().neighbor(grid().tree().neighbor(cellId, face[1]), face[0]);
5621 reverseEdge[2] = reverseEdgeStencil[1][edge];
5622 reverseEdge[3] = reverseEdgeStencil[2][edge];
5629 const MBool isGapCell = candidates[cndt].isGapCell;
5630 if(m_complexBoundary && m_noLevelSetsUsedForMb > 1 && (!isGapCell)) {
5632 endSet = m_noLevelSetsUsedForMb;
5636 for(
MInt set = startSet; set < endSet; set++) {
5639 phi[0] = candidates[cndt].nodalValues[set][nodeStencil[0][edge]];
5640 phi[1] = candidates[cndt].nodalValues[set][nodeStencil[1][edge]];
5642 ASSERT(candidates[cndt].nodeValueSet[nodeStencil[0][edge]],
"");
5643 ASSERT(candidates[cndt].nodeValueSet[nodeStencil[1][edge]],
"");
5646 if(fabs(phi[0]) < eps) {
5652 if(fabs(phi[1]) < eps) {
5662 if(m_multiCutCell) {
5663 if(
approx(fabs(phi[0]), cellHalfLength, eps)) {
5665 phi[0] = cellHalfLength - eps;
5667 phi[0] = -cellHalfLength + eps;
5669 if(
approx(fabs(phi[1]), cellHalfLength, eps)) {
5671 phi[1] = cellHalfLength - eps;
5673 phi[1] = -cellHalfLength + eps;
5678 if(phi[0] * phi[1] < 0) {
5681 const MFloat deltaCutpoint = (phi[0] + phi[1]) / (phi[0] - phi[1]);
5682 ASSERT(deltaCutpoint < F1 && deltaCutpoint > -F1,
"");
5683 for(
MInt dim = 0; dim < nDim_; dim++) {
5684 if(dim == edgeDOF)
continue;
5685 cutpoint[dim] = a_coordinate(cellId, dim) + signStencil[dim][edge] * cellHalfLength;
5687 cutpoint[edgeDOF] = a_coordinate(cellId, edgeDOF) + cellHalfLength * deltaCutpoint;
5690 for(
MInt nb = 0; nb < edgesPerDim; nb++) {
5691 const MInt nbCell = directNeighbourIds[nb];
5692 if(nbCell < 0)
continue;
5693 const MInt nbCandidate = candidateIds[nbCell];
5697 if(nbCandidate < 0)
continue;
5700 ASSERT(!candidates[nbCandidate].edgeChecked[reverseEdge[nb]],
"");
5704 candidatesOrder.push_back(nbCandidate);
5707 const MInt noCPs = candidates[nbCandidate].noCutPoints;
5708 if(noCPs == 2 * m_noEdges) {
5709 mTerm(1, AT_,
"Error: too many cut points, can't store more!");
5711 if(candidates[nbCandidate].noCutPointsOnEdge[reverseEdge[nb]] == 2) {
5713 mTerm(1, AT_,
"Error: already two cut points on edge, can't store more!");
5722 for(
MInt dim = 0; dim < nDim_; dim++) {
5723 candidates[nbCandidate].cutPoints[noCPs][dim] = cutpoint[dim];
5727 ASSERT(candidates[nbCandidate].associatedBodyIds[set] > -1,
"cutPoint with invalid bodyId! ");
5729 candidates[nbCandidate].cutBodyIds[noCPs] = candidates[nbCandidate].associatedBodyIds[set];
5730 candidates[nbCandidate].cutEdges[noCPs] = reverseEdge[nb];
5731 candidates[nbCandidate].noCutPoints++;
5732 candidates[nbCandidate].noCutPointsOnEdge[reverseEdge[nb]]++;
5739 for(
MInt nb = 0; nb < edgesPerDim; nb++) {
5740 if(directNeighbourIds[nb] < 0)
continue;
5741 const MInt nbCndt = candidateIds[directNeighbourIds[nb]];
5742 if(nbCndt < 0)
continue;
5743 candidates[nbCndt].edgeChecked[reverseEdge[nb]] =
true;
5749 MPI_Allreduce(MPI_IN_PLACE, &errorFlag, 1, MPI_INT, MPI_MAX, grid().mpiComm(), AT_,
"MPI_IN_PLACE",
"errorFlag");
5752 mTerm(1, AT_,
"Critical error in computeCutPoints! Quit.");
5767template <MInt nDim_>
5769 const MFloat* scalarField,
const MInt* bodyIdField,
5770 const MBool*
const gapPropertyField,
const MBool gapClosure) {
5774 const MInt debugTimeStep = -2;
5775 const MInt debugGlobalId = 216221;
5778 m_cutLvlJumpCandidates.clear();
5780 if(candidates.empty())
return;
5782 const MInt noCandidates = (signed)candidates.size();
5785 if(m_bndryLvlJumps) {
5786 vector<std::pair<MInt, MInt>> diagonalJumpCells;
5787 for(
MInt cnd = 0; cnd < noCandidates; cnd++) {
5788 const MInt cellId = candidates[cnd].cellId;
5790 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
5791 cerr <<
"Cell is candidate!" << endl;
5794 ASSERT(!candidates[cnd].isbndryLvlJumpParent,
"");
5796 diagonalJumpCells.clear();
5797 const MInt parentId = grid().tree().parent(cellId);
5798 if(parentId < 0)
continue;
5799 for(
MInt dir = 0; dir < m_noDirs; dir++) {
5801 if(!grid().tree().hasNeighbor(cellId, dir) && grid().tree().hasNeighbor(parentId, dir)) {
5804 id = m_cutLvlJumpCandidates.size();
5805 m_cutLvlJumpCandidates.emplace_back();
5806 m_cutLvlJumpCandidates[
id].candId = cnd;
5809 for(child = 0; child < m_maxNoChilds; child++) {
5810 const MInt childCellId = grid().tree().child(parentId, child);
5811 if(childCellId == cellId)
break;
5813 m_cutLvlJumpCandidates[
id].childId = child;
5815 if(candidateIds[parentId] < 0) {
5816 const MInt newId = candidates.size();
5817 candidates.emplace_back();
5818 candidates[newId].cellId = parentId;
5819 candidates[newId].isbndryLvlJumpParent =
true;
5820 candidateIds[parentId] = newId;
5822 if(!candidates[candidateIds[parentId]].isbndryLvlJumpParent) {
5823 candidates[candidateIds[parentId]].isbndryLvlJumpParent =
true;
5826 const MInt parentCandId = candidateIds[parentId];
5827 ASSERT(parentCandId > -1,
"");
5828 m_cutLvlJumpCandidates[
id].parentCandId = parentCandId;
5831 const MInt noJumps = m_cutLvlJumpCandidates[
id].noJumps;
5834 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
5835 cerr <<
"Direct nieghbor Jump!" << dir << endl;
5839 m_cutLvlJumpCandidates[
id].dirs[noJumps] = dir;
5840 m_cutLvlJumpCandidates[
id].diagonalDirs[noJumps] = -2;
5841 m_cutLvlJumpCandidates[
id].neighborType[noJumps] = 0;
5843 m_cutLvlJumpCandidates[
id].noJumps++;
5846 const MInt nghbrId1 = grid().tree().neighbor(cellId, dir);
5847 if(nghbrId1 < 0)
continue;
5848 if(!grid().tree().isLeafCell(nghbrId1))
continue;
5849 const MInt ngbParent1 = grid().tree().parent(nghbrId1);
5850 if(ngbParent1 < 0)
continue;
5851 if(ngbParent1 == parentId)
continue;
5852 for(
MInt dir1 = 0; dir1 < m_noDirs; dir1++) {
5853 if((dir1 / 2) == (dir / 2))
continue;
5855 if(!grid().tree().hasNeighbor(nghbrId1, dir1) && grid().tree().hasNeighbor(ngbParent1, dir1)) {
5857 MBool alreadyAdded =
false;
5859 for(i = 0; i < diagonalJumpCells.size(); i++) {
5860 if(diagonalJumpCells[i].first == grid().tree().neighbor(ngbParent1, dir1)
5861 && diagonalJumpCells[i].second == 1) {
5862 alreadyAdded =
true;
5864 }
else if(diagonalJumpCells[i].first == grid().tree().neighbor(ngbParent1, dir1)) {
5871 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
5872 cerr <<
"Already found: " << dir <<
" " << dir1 <<
" "
5873 << grid().tree().globalId(grid().tree().neighbor(ngbParent1, dir1)) << endl;
5877 ASSERT(
id > -1,
"");
5882 diagonalJumpCells.emplace_back(make_pair(grid().tree().neighbor(ngbParent1, dir1), 1));
5886 id = m_cutLvlJumpCandidates.size();
5887 m_cutLvlJumpCandidates.emplace_back();
5888 m_cutLvlJumpCandidates[
id].candId = cnd;
5891 for(child = 0; child < m_maxNoChilds; child++) {
5892 const MInt childCellId = grid().tree().child(parentId, child);
5893 if(childCellId == cellId)
break;
5895 m_cutLvlJumpCandidates[
id].childId = child;
5897 if(candidateIds[parentId] < 0) {
5898 const MInt newId = candidates.size();
5899 candidates.emplace_back();
5900 candidates[newId].cellId = parentId;
5901 candidates[newId].isbndryLvlJumpParent =
true;
5902 candidateIds[parentId] = newId;
5904 if(!candidates[candidateIds[parentId]].isbndryLvlJumpParent) {
5905 candidates[candidateIds[parentId]].isbndryLvlJumpParent =
true;
5908 const MInt parentCandId = candidateIds[parentId];
5909 ASSERT(parentCandId > -1,
"");
5910 m_cutLvlJumpCandidates[
id].parentCandId = parentCandId;
5912 const MInt noJumps = m_cutLvlJumpCandidates[
id].noJumps;
5915 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
5916 cerr <<
"Adding type1: " << dir <<
" " << dir1 <<
" "
5917 << grid().tree().globalId(grid().tree().neighbor(ngbParent1, dir1)) <<
" " << noJumps <<
" "
5922 m_cutLvlJumpCandidates[
id].neighborType[noJumps] = 1;
5924 m_cutLvlJumpCandidates[
id].dirs[noJumps] = dir;
5925 m_cutLvlJumpCandidates[
id].diagonalDirs[noJumps] = dir1;
5927 m_cutLvlJumpCandidates[
id].noJumps++;
5930 IF_CONSTEXPR(nDim == 3) {
5931 const MInt nghbrId2 = grid().tree().neighbor(nghbrId1, dir1);
5932 if(nghbrId2 < 0)
continue;
5933 if(!grid().tree().isLeafCell(nghbrId2))
continue;
5934 const MInt ngbParent2 = grid().tree().parent(nghbrId2);
5935 if(ngbParent2 < 0)
continue;
5936 if(ngbParent2 == parentId)
continue;
5937 for(
MInt dir2 = 0; dir2 < 2 * nDim; dir2++) {
5938 if(((dir2 / 2) == (dir / 2)) || ((dir2 / 2) == (dir1 / 2)))
continue;
5940 if(!grid().tree().hasNeighbor(nghbrId2, dir2) && grid().tree().hasNeighbor(ngbParent2, dir2)) {
5942 MBool alreadyAdded =
false;
5943 for(
MUint i = 0; i < diagonalJumpCells.size(); i++) {
5944 if(diagonalJumpCells[i].first == grid().tree().neighbor(ngbParent2, dir2)) {
5945 alreadyAdded =
true;
5950 ASSERT(
id > -1,
"");
5954 diagonalJumpCells.emplace_back(make_pair(grid().tree().neighbor(ngbParent2, dir2), 2));
5957 id = m_cutLvlJumpCandidates.size();
5958 m_cutLvlJumpCandidates.emplace_back();
5959 m_cutLvlJumpCandidates[
id].candId = cnd;
5962 for(child = 0; child < m_maxNoChilds; child++) {
5963 const MInt childCellId = grid().tree().child(parentId, child);
5964 if(childCellId == cellId)
break;
5966 m_cutLvlJumpCandidates[
id].childId = child;
5968 if(candidateIds[parentId] < 0) {
5969 const MInt newId = candidates.size();
5970 candidates.emplace_back();
5971 candidates[newId].cellId = parentId;
5972 candidates[newId].isbndryLvlJumpParent =
true;
5973 candidateIds[parentId] = newId;
5975 if(!candidates[candidateIds[parentId]].isbndryLvlJumpParent) {
5976 candidates[candidateIds[parentId]].isbndryLvlJumpParent =
true;
5979 const MInt parentCandId = candidateIds[parentId];
5980 ASSERT(parentCandId > -1,
"");
5981 m_cutLvlJumpCandidates[
id].parentCandId = parentCandId;
5983 const MInt noJumps = m_cutLvlJumpCandidates[
id].noJumps;
5984 m_cutLvlJumpCandidates[
id].neighborType[noJumps] = 2;
5986 m_cutLvlJumpCandidates[
id].dirs[noJumps] = dir;
5987 m_cutLvlJumpCandidates[
id].diagonalDirs[noJumps] = dir1;
5988 m_cutLvlJumpCandidates[
id].noJumps++;
5991 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
5992 cerr <<
"Adding type2: " << dir <<
" " << dir1 <<
" " << dir2 <<
" "
5993 << grid().tree().globalId(grid().tree().neighbor(ngbParent2, dir1)) <<
" " << noJumps
6006 ASSERT(noCandidates <= (
signed)candidates.size(),
"");
6008 const MInt nodeStencil[3][8] = {{0, 1, 0, 1, 0, 1, 0, 1}, {2, 2, 3, 3, 2, 2, 3, 3}, {4, 4, 4, 4, 5, 5, 5, 5}};
6010 const MInt reverseNode[3][8] = {{1, 0, 3, 2, 5, 4, 7, 6}, {2, 3, 0, 1, 6, 7, 4, 5}, {4, 5, 6, 7, 0, 1, 2, 3}};
6012 MInt nghbrIds[m_noCorners]{};
6013 MInt nghbrNodes[m_noCorners]{};
6015 for(
MInt cnd = 0; cnd < (signed)candidates.size(); cnd++) {
6016 const MInt cellId = candidates[cnd].cellId;
6019 candidates[cnd].isGapCell = gapPropertyField[cellId];
6020 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6021 if(bodyIdField !=
nullptr) {
6022 candidates[cnd].associatedBodyIds[set] = bodyIdField[IDX_LSSETMB(cellId, set)];
6026 for(
MInt node = 0; node < m_noCorners; node++) {
6027 if(candidates[cnd].nodeValueSet[node])
continue;
6030 MInt noNeighborsPerNode = 0;
6032 nghbrIds[noNeighborsPerNode] = cellId;
6033 nghbrNodes[noNeighborsPerNode] = node;
6034 noNeighborsPerNode++;
6036 for(
MInt i = 0; i < nDim; i++) {
6037 const MInt firstDir = nodeStencil[i][node];
6038 const MInt firstNghbrNode = reverseNode[i][node];
6039 const MInt firstNghbrId = grid().tree().neighbor(cellId, firstDir);
6040 if(firstNghbrId > -1) {
6041 nghbrIds[noNeighborsPerNode] = firstNghbrId;
6042 nghbrNodes[noNeighborsPerNode] = firstNghbrNode;
6043 noNeighborsPerNode++;
6044 for(
MInt j = i + 1; j < nDim; j++) {
6045 const MInt secondDir = nodeStencil[j][node];
6046 if(secondDir == firstDir)
continue;
6047 const MInt secondNghbrNode = reverseNode[j][firstNghbrNode];
6048 const MInt secondNghbrId = grid().tree().neighbor(firstNghbrId, secondDir);
6049 if(secondNghbrId > -1) {
6050 nghbrIds[noNeighborsPerNode] = secondNghbrId;
6051 nghbrNodes[noNeighborsPerNode] = secondNghbrNode;
6052 noNeighborsPerNode++;
6053 IF_CONSTEXPR(nDim == 3) {
6054 for(
MInt k = j + 1; k < nDim; k++) {
6055 const MInt thirdDir = nodeStencil[k][node];
6056 if(thirdDir == firstDir || thirdDir == secondDir)
continue;
6057 const MInt thirdNghbrNode = reverseNode[k][secondNghbrNode];
6058 const MInt thirdNghbrId = grid().tree().neighbor(secondNghbrId, thirdDir);
6059 if(thirdNghbrId > -1) {
6060 nghbrIds[noNeighborsPerNode] = thirdNghbrId;
6061 nghbrNodes[noNeighborsPerNode] = thirdNghbrNode;
6062 noNeighborsPerNode++;
6071 MBool gapBoundary =
false;
6072 MBool setBoundary =
false;
6080 const MBool isGapCell = gapPropertyField[cellId];
6081 for(
MInt nghbrNode = 0; nghbrNode < noNeighborsPerNode; nghbrNode++) {
6082 const MBool isGapNghbr = gapPropertyField[nghbrIds[nghbrNode]];
6084 if(isGapCell != isGapNghbr) {
6088 const MInt bodyId = bodyIdField[IDX_LSSETMB(cellId, 0)];
6089 for(
MInt nghbrNode = 0; nghbrNode < noNeighborsPerNode; nghbrNode++) {
6090 const MInt bodyIdNghbr = bodyIdField[IDX_LSSETMB(nghbrIds[nghbrNode], 0)];
6091 if(bodyId != bodyIdNghbr) setBoundary =
true;
6095 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6098 for(
MInt nghbrNode = 0; nghbrNode < noNeighborsPerNode; nghbrNode++) {
6102 const MInt bodyIdNghbr = bodyIdField[IDX_LSSETMB(nghbrIds[nghbrNode], 0)];
6106 if(m_bodyToSetTable[bodyIdNghbr] == set) usedSet = 0;
6116 phi += scalarField[IDX_LSSETMB(nghbrIds[nghbrNode], usedSet)];
6119 phi /= noNeighborsPerNode;
6121 for(
MInt nghbrNode = 0; nghbrNode < noNeighborsPerNode; nghbrNode++) {
6122 const MInt nghbrCand = candidateIds[nghbrIds[nghbrNode]];
6123 if(nghbrCand < 0)
continue;
6124 candidates[nghbrCand].nodeValueSet[nghbrNodes[nghbrNode]] =
true;
6125 candidates[nghbrCand].nodalValues[set][nghbrNodes[nghbrNode]] = phi;
6130 for(
MInt nghbrNode = 0; nghbrNode < noNeighborsPerNode; nghbrNode++) {
6132 const MInt candidateId = candidateIds[nghbrIds[nghbrNode]];
6133 if(candidateId < 0)
continue;
6134 const MInt cellId2 = candidates[candidateId].cellId;
6135 ASSERT(cellId2 == nghbrIds[nghbrNode],
"");
6136 const MInt set = m_bodyToSetTable[bodyIdField[IDX_LSSETMB(cellId2, 0)]];
6137 if(abs(candidates[candidateId].nodalValues[0][nghbrNodes[nghbrNode]]
6138 - candidates[candidateId].nodalValues[set][nghbrNodes[nghbrNode]])
6140 cerr <<
"ERROR in Node-LS-Values " << cellId2 <<
" globalId " << grid().tree().globalId(cellId2) <<
" "
6141 << grid().tree().globalId(cellId) <<
" "
6143 <<
" Gap " << candidates[cnd].isGapCell <<
" " << candidates[candidateId].isGapCell <<
" BodyId "
6144 << candidates[cnd].associatedBodyIds[0] <<
" " << candidates[candidateId].associatedBodyIds[0] <<
" "
6145 << nghbrNode << endl;
6152 if(m_bndryLvlJumps && !m_cutLvlJumpCandidates.empty()) {
6153 correctNodalValuesAtLevelJump(candidates, candidateIds);
6157 for(
MInt cnd = noCandidates; cnd < (signed)candidates.size(); cnd++) {
6158 const MInt parentId = candidates[cnd].cellId;
6159 ASSERT(candidates[cnd].isbndryLvlJumpParent,
"");
6160 candidateIds[parentId] = -1;
6163 candidates.resize(noCandidates);
6171template <MInt nDim_>
6173 const MInt* noMaxLevelWindowCells,
6174 const MInt** maxLevelHaloCells,
6176 MInt* candidateIds) {
6179 auto* mpi_send_req =
new MPI_Request[grid().noNeighborDomains()];
6180 auto* mpi_recv_req =
new MPI_Request[grid().noNeighborDomains()];
6182 MIntScratchSpace sendBufferSize(grid().noNeighborDomains(), AT_,
"sendBufferSize");
6183 MIntScratchSpace receiveBufferSize(grid().noNeighborDomains(), AT_,
"receiveBufferSize");
6185 MIntScratchSpace sendBuffersInt(grid().noNeighborDomains(), AT_,
"sendBuffersInt");
6186 MIntScratchSpace receiveBuffersInt(grid().noNeighborDomains(), AT_,
"receiveBuffersInt");
6188 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6189 mpi_send_req[i] = MPI_REQUEST_NULL;
6190 mpi_recv_req[i] = MPI_REQUEST_NULL;
6194 MInt sendBufferCounter;
6195 MInt sendBufferCounterOverall = 0;
6196 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6197 sendBufferCounter = 0;
6198 sendBufferCounter++;
6199 for(
MInt j = 0; j < noMaxLevelWindowCells[i]; j++) {
6200 const MInt cellId = maxLevelWindowCells[i][j];
6201 const MInt cndId = candidateIds[cellId];
6202 if(cndId < 0)
continue;
6203 for(
MInt node = 0; node < m_noCorners; node++) {
6204 ASSERT(candidates[cndId].nodeValueSet[node],
"");
6208 sendBufferCounter++;
6211 sendBufferCounter += m_noLevelSetsUsedForMb * m_noCorners;
6214 sendBufferCounter++;
6217 sendBufferCounter += m_noLevelSetsUsedForMb;
6219 sendBufferSize[i] = sendBufferCounter;
6220 sendBuffersInt[i] = sendBufferCounter;
6221 sendBufferCounterOverall += sendBufferCounter;
6225 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6226 MPI_Irecv(&receiveBuffersInt[i], 1, MPI_INT, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_recv_req[i], AT_,
6227 "receiveBuffers[i]");
6229 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6230 MPI_Isend(&sendBuffersInt[i], 1, MPI_INT, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_send_req[i], AT_,
6231 "sendBuffersInt[i]");
6233 MPI_Waitall(grid().noNeighborDomains(), mpi_recv_req, MPI_STATUSES_IGNORE, AT_);
6234 MPI_Waitall(grid().noNeighborDomains(), mpi_send_req, MPI_STATUSES_IGNORE, AT_);
6237 MInt receiveBufferCounterOverall = 0;
6238 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6239 receiveBufferCounterOverall += receiveBuffersInt[i];
6240 receiveBufferSize[i] = receiveBuffersInt[i];
6243 MFloatScratchSpace sendBuffersOverall(sendBufferCounterOverall, AT_,
"sendBuffersOverall");
6244 MFloatScratchSpace receiveBuffersOverall(receiveBufferCounterOverall, AT_,
"receiveBuffersOverall");
6247 MInt counterSend = 0;
6248 MInt counterReceive = 0;
6249 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6250 sendBuffers.
p[i] = &sendBuffersOverall.
p[counterSend];
6251 counterSend += sendBufferSize[i];
6252 receiveBuffers.
p[i] = &receiveBuffersOverall.
p[counterReceive];
6253 counterReceive += receiveBufferSize[i];
6257 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6258 mpi_send_req[i] = MPI_REQUEST_NULL;
6259 mpi_recv_req[i] = MPI_REQUEST_NULL;
6262 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6263 sendBufferCounter = 0;
6264 sendBuffers[i][0] = (
MFloat)(-1);
6265 sendBufferCounter++;
6266 for(
MInt j = 0; j < noMaxLevelWindowCells[i]; j++) {
6267 const MInt cellId = maxLevelWindowCells[i][j];
6268 const MInt cndId = candidateIds[cellId];
6269 if(cndId < 0)
continue;
6271 sendBuffers[i][sendBufferCounter] = (
MFloat)j;
6272 sendBufferCounter++;
6274 for(
MInt node = 0; node < m_noCorners; node++) {
6275 ASSERT(candidates[cndId].nodeValueSet[node],
"");
6276 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6277 sendBuffers[i][sendBufferCounter] = candidates[cndId].nodalValues[set][node];
6278 sendBufferCounter++;
6282 sendBuffers[i][sendBufferCounter] = (
MFloat)candidates[cndId].isGapCell;
6283 sendBufferCounter++;
6285 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6286 sendBuffers[i][sendBufferCounter] = (
MFloat)candidates[cndId].associatedBodyIds[set];
6287 sendBufferCounter++;
6290 sendBufferSize[i] = sendBufferCounter;
6291 sendBuffers[i][0] = (
MFloat)(sendBufferCounter);
6296 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6297 MInt bufSize = receiveBufferSize[i];
6298 MPI_Irecv(receiveBuffers[i], bufSize, MPI_DOUBLE, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_recv_req[i],
6299 AT_,
"receiveBuffers[i]");
6303 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6304 MInt bufSize = sendBufferSize[i];
6305 MPI_Isend(sendBuffers[i], bufSize, MPI_DOUBLE, grid().neighborDomain(i), 2, grid().mpiComm(), &mpi_send_req[i], AT_,
6308 MPI_Waitall(grid().noNeighborDomains(), mpi_recv_req, MPI_STATUSES_IGNORE, AT_);
6309 MPI_Waitall(grid().noNeighborDomains(), mpi_send_req, MPI_STATUSES_IGNORE, AT_);
6313 MInt receiveBufferCounter = 0;
6315 for(
MInt i = 0; i < grid().noNeighborDomains(); i++) {
6316 receiveBufferCounter = 0;
6317 if(receiveBufferSize[i] != receiveBuffersInt[i]) {
6318 mTerm(1, AT_,
" receiveBufferSize doesn't match expected size from previous communication! ");
6320 if(receiveBufferSize[i] == 0) {
6321 m_log <<
"Warning: empty message from rank " << grid().neighborDomain(i) << endl;
6323 receiveBufferCounter++;
6324 while(receiveBufferCounter < receiveBufferSize[i]) {
6325 j = (
MInt)receiveBuffers[i][receiveBufferCounter];
6326 receiveBufferCounter++;
6327 const MInt cellId = maxLevelHaloCells[i][j];
6330 if(cellId > grid().tree().size()) skip =
true;
6333 const MInt candId = candidates.size();
6334 candidates.emplace_back();
6335 candidates[candId].cellId = cellId;
6338 for(
MInt node = 0; node < m_noCorners; node++) {
6339 candidates[candId].nodeValueSet[node] =
true;
6340 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6341 candidates[candId].nodalValues[set][node] = receiveBuffers[i][receiveBufferCounter];
6342 receiveBufferCounter++;
6346 candidates[candId].isGapCell = (
MInt)receiveBuffers[i][receiveBufferCounter];
6347 receiveBufferCounter++;
6349 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6350 candidates[candId].associatedBodyIds[set] = (
MInt)receiveBuffers[i][receiveBufferCounter];
6351 receiveBufferCounter++;
6354 candidateIds[cellId] = candId;
6358 receiveBufferCounter += m_noLevelSetsUsedForMb * m_noCorners;
6363 delete[] mpi_send_req;
6364 delete[] mpi_recv_req;
6380 const MInt faceCornerMapping[6][4] = {{2, 6, 4, 0}, {1, 5, 7, 3}, {0, 4, 5, 1},
6381 {3, 7, 6, 2}, {2, 3, 1, 0}, {6, 7, 5, 4}};
6384 ASSERT(m_noLevelSetsUsedForMb <= CC::maxNoSets,
"");
6386 for(
MInt cndt = 0; cndt < (signed)candidates.size(); cndt++) {
6387 if(candidates[cndt].noCutPoints < nDim)
continue;
6389 const MInt cutCell = cutCellData.size();
6390 const MInt cellId = candidates[cndt].cellId;
6392 cutCellIdMapping[cellId] = cndt;
6394 cutCellData.emplace_back();
6395 cutCellData[cutCell].cellId = cellId;
6396 cutCellData[cutCell].isGapCell = candidates[cndt].isGapCell;
6397 ASSERT(candidates[cndt].noCutPoints <= CC::maxNoCutPoints,
"");
6398 cutCellData[cutCell].noCutPoints = candidates[cndt].noCutPoints;
6399 for(
MInt cp = 0; cp < candidates[cndt].noCutPoints; cp++) {
6400 for(
MInt i = 0; i < nDim; i++) {
6401 cutCellData[cutCell].cutPoints[cp][i] = candidates[cndt].cutPoints[cp][i];
6403 cutCellData[cutCell].cutBodyIds[cp] = candidates[cndt].cutBodyIds[cp];
6404 cutCellData[cutCell].cutEdges[cp] = candidates[cndt].cutEdges[cp];
6407 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6408 cutCellData[cutCell].associatedBodyIds[set] = candidates[cndt].associatedBodyIds[set];
6410 for(
MInt k = 0; k < CC::noFaces; k++) {
6412 for(
MInt j = 0; j < 4; j++) {
6413 phi += F1B4 * candidates[cndt].nodalValues[set][faceCornerMapping[k][j]];
6415 cutCellData[cutCell].faceCentroidIsInsideGeometry[set][k] = (phi < F0);
6421 if(m_complexBoundary && m_noLevelSetsUsedForMb > 1 && (!cutCellData[cutCell].isGapCell)) {
6423 endSet = m_noLevelSetsUsedForMb;
6426 for(
MInt node = 0; node < m_noCorners; node++) {
6427 ASSERT(candidates[cndt].nodeValueSet[node],
"");
6428 MBool isInside =
false;
6429 for(
MInt set = startSet; set < endSet; set++) {
6430 const MBool isInsideSet = candidates[cndt].nodalValues[set][node] < F0;
6431 isInside = isInside || isInsideSet;
6432 cutCellData[cutCell].cornerIsInsideGeometry[set][node] = isInsideSet;
6434 cutCellData[cutCell].cornerIsInsideGeometry[0][node] = isInside;
6447template <MInt nDim_>
6451 if(candidates.empty())
return;
6454 m_scaledCutCell =
false;
6456 const MInt DOFStencil[12] = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
6458 const MFloat signStencil[8][3] = {{-F1, -F1, -F1}, {F1, -F1, -F1}, {-F1, F1, -F1}, {F1, F1, -F1},
6459 {-F1, -F1, F1}, {F1, -F1, F1}, {-F1, F1, F1}, {F1, F1, F1}};
6461 const MInt nodeStencil[2][12] = {{0, 1, 0, 2, 4, 5, 4, 6, 0, 1, 2, 3}, {2, 3, 1, 3, 6, 7, 5, 7, 4, 5, 6, 7}};
6463 const MFloat epsilon = 0.00000000001;
6464 const MInt maxNoCutPointsPerEdge = 10;
6465 std::array<MFloat, nDim>
a;
6466 std::array<MFloat, nDim>
b;
6467 std::array<MFloat, nDim> c;
6468 std::array<MFloat, nDim> d;
6469 std::array<MFloat, nDim> e;
6470 std::array<MFloat, nDim> pP;
6473 MFloat corners[m_noCorners][nDim];
6479 const MInt m_maxNoBndryCndIds = 100;
6481 MInt noBndryCndIds = 0;
6482 MInt bndryCndIds[m_maxNoBndryCndIds]{};
6486 for(
MInt cndt = 0; cndt < (signed)candidates.size(); cndt++) {
6487 ASSERT(candidates[cndt].noCutPoints == 0,
"");
6488 const MInt cellId = candidates[cndt].cellId;
6490 candidates[cndt].associatedBodyIds[0] = std::numeric_limits<MInt>::max();
6492 ASSERT(cellId < grid().tree().size(),
"");
6494 if(grid().tree().noChildren(cellId) > 0)
continue;
6495 const MFloat cellLength = a_cellLengthAtCell(cellId);
6496 const MFloat cellHalfLength = F1B2 * cellLength;
6498 const MFloat eps = cellLength / 10000000.0;
6501 for(
MInt node = 0; node < m_noCorners; node++) {
6502 for(
MInt i = 0; i < nDim; i++) {
6503 corners[node][i] = grid().tree().coordinate(cellId, i) + signStencil[node][i] * cellHalfLength;
6508 for(
MInt i = 0; i < nDim; i++) {
6509 target[i] = grid().tree().coordinate(cellId, i) - cellHalfLength * 1.05;
6510 target[i + nDim] = grid().tree().coordinate(cellId, i) + cellHalfLength * 1.05;
6514 std::vector<MInt> nodeList;
6515 if(m_gridCutTest ==
"SAT") {
6516 geometry().getIntersectionElements(target, nodeList, cellHalfLength, &a_coordinate(cellId, 0));
6518 geometry().getIntersectionElements(target, nodeList);
6523 for(
MInt edge = 0; edge < m_noEdges; edge++) {
6524 const MInt edgeDOF = DOFStencil[edge];
6526 ASSERT(candidates[cndt].noCutPointsOnEdge[edge] == 0,
"");
6528 MBool edgeCut =
false;
6531 for(
MInt i = 0; i < nDim; i++) {
6532 d[i] = corners[nodeStencil[0][edge]][i];
6533 e[i] = corners[nodeStencil[1][edge]][i];
6537 for(
MInt n = 0; n < (signed)nodeList.size(); n++) {
6538 MBool cutsEdge =
false;
6540 IF_CONSTEXPR(nDim == 3) {
6541 const MInt spaceId = (edgeDOF + 1) % nDim;
6542 const MInt spaceId1 = (edgeDOF + 2) % nDim;
6543 const MInt spaceId2 = edgeDOF;
6547 for(
MInt k = 0; k < nDim; k++) {
6548 a[k] = m_geometry->elements[nodeList[n]].m_vertices[0][k];
6549 b[k] = m_geometry->elements[nodeList[n]].m_vertices[1][k];
6550 c[k] = m_geometry->elements[nodeList[n]].m_vertices[2][k];
6552 if(
approx(
a[spaceId1],
b[spaceId1], MFloatEps) && !
approx(
a[spaceId1], c[spaceId1], MFloatEps)) {
6554 q = (d[spaceId1] -
a[spaceId1]) / (c[spaceId1] -
a[spaceId1]);
6555 p = (d[spaceId] -
a[spaceId] - q * (c[spaceId] -
a[spaceId])) / (
b[spaceId] -
a[spaceId]);
6557 if(!
approx(
a[spaceId1],
b[spaceId1], MFloatEps) &&
approx(
a[spaceId1], c[spaceId1], MFloatEps)) {
6559 p = (d[spaceId1] -
a[spaceId1]) / (
b[spaceId1] -
a[spaceId1]);
6560 q = (d[spaceId] -
a[spaceId] - p * (
b[spaceId] -
a[spaceId])) / (c[spaceId] -
a[spaceId]);
6562 if(
approx(
a[spaceId],
b[spaceId], MFloatEps) && !
approx(
a[spaceId], c[spaceId], MFloatEps)) {
6564 q = (d[spaceId] -
a[spaceId]) / (c[spaceId] -
a[spaceId]);
6565 p = (d[spaceId1] -
a[spaceId1] - q * (c[spaceId1] -
a[spaceId1])) / (
b[spaceId1] -
a[spaceId1]);
6567 if(!
approx(
a[spaceId],
b[spaceId], MFloatEps) &&
approx(
a[spaceId], c[spaceId], MFloatEps)) {
6569 p = (d[spaceId] -
a[spaceId]) / (
b[spaceId] -
a[spaceId]);
6570 q = (d[spaceId1] -
a[spaceId1] - p * (
b[spaceId1] -
a[spaceId1])) / (c[spaceId1] -
a[spaceId1]);
6573 q = ((d[spaceId1] -
a[spaceId1]) * (
b[spaceId] -
a[spaceId])
6574 - (
b[spaceId1] -
a[spaceId1]) * (d[spaceId] -
a[spaceId]))
6575 / ((c[spaceId1] -
a[spaceId1]) * (
b[spaceId] -
a[spaceId])
6576 - (
b[spaceId1] -
a[spaceId1]) * (c[spaceId] -
a[spaceId]));
6577 p = (d[spaceId] -
a[spaceId] - q * (c[spaceId] -
a[spaceId])) / (
b[spaceId] -
a[spaceId]);
6584 if(p * q >= 0 || p * q < 0) {
6586 MFloat gamma =
a[spaceId2] + p * (
b[spaceId2] -
a[spaceId2]) + q * (c[spaceId2] -
a[spaceId2]);
6587 MFloat s = (gamma - d[spaceId2]) / (e[spaceId2] - d[spaceId2]);
6589 if(s < -epsilon || s > F1 + epsilon || p < -epsilon || q < -epsilon || (p + q) > F1 + epsilon) {
6593 if(candidates[cndt].noCutPointsOnEdge[edge] < maxNoCutPointsPerEdge) {
6594 for(
MInt k = 0; k < nDim; k++) {
6595 pP[k] = d[k] + s * (e[k] - d[k]);
6596 cutPointsEdge(candidates[cndt].noCutPointsOnEdge[edge], k) = pP[k];
6599 mTerm(1, AT_,
" Too many cut points on edge...");
6605 if(geometry().edgeTriangleIntersection(geometry().elements[nodeList[n]].m_vertices[0],
6606 geometry().elements[nodeList[n]].m_vertices[1],
nullptr, d.data(),
6608 for(
MInt k = 0; k < nDim; k++) {
6609 a[k] = geometry().elements[nodeList[n]].m_vertices[0][k];
6610 b[k] = geometry().elements[nodeList[n]].m_vertices[1][k];
6613 MFloat gamma = (
b[0] -
a[0]) * (d[1] - e[1]) - (d[0] - e[0]) * (
b[1] -
a[1]);
6614 if(
ABS(gamma) < 0.0000000000001)
continue;
6615 MFloat s1 = ((d[0] - e[0]) * (
a[1] - d[1]) - (d[1] - e[1]) * (
a[0] - d[0])) / gamma;
6616 MFloat s2 = ((
b[0] -
a[0]) * (d[1] -
a[1]) - (d[0] -
a[0]) * (
b[1] -
a[1])) / gamma;
6621 if(candidates[cndt].noCutPointsOnEdge[edge] < maxNoCutPointsPerEdge) {
6622 for(
MInt k = 0; k < nDim; k++) {
6623 if(s1 * s1 < s2 * s2)
6624 pP[k] = d[k] + s2 * (e[k] - d[k]);
6626 pP[k] =
a[k] + s1 * (
b[k] -
a[k]);
6627 cutPointsEdge(candidates[cndt].noCutPointsOnEdge[edge], k) = pP[k];
6630 mTerm(1, AT_,
" Too many cut points on edge...");
6637 if(candidates[cndt].noCutPoints < m_noEdges) {
6639 const MInt bndCndId = geometry().elements[nodeList[n]].m_bndCndId;
6640 if(bndCndId < candidates[cndt].associatedBodyIds[0]) {
6641 candidates[cndt].associatedBodyIds[0] = bndCndId;
6642 MBool newBndryCnd =
true;
6643 for(
MInt j = 0; j < noBndryCndIds; j++) {
6644 if(bndryCndIds[j] == bndCndId) {
6645 newBndryCnd =
false;
6650 ASSERT(noBndryCndIds < m_maxNoBndryCndIds,
"Too many different boundary conditions...");
6651 bndryCndIds[noBndryCndIds] = bndCndId;
6659 MBool newCutPoint =
true;
6660 for(
MInt i = 0; i < candidates[cndt].noCutPointsOnEdge[edge]; i++) {
6661 MBool equal = (
ABS(pP[0] - cutPointsEdge(i, 0)) < eps &&
ABS(pP[1] - cutPointsEdge(i, 1)) < eps);
6662 IF_CONSTEXPR(nDim == 3) equal = (equal && (
ABS(pP[2] - cutPointsEdge(i, 2)) < eps));
6664 newCutPoint =
false;
6670 candidates[cndt].noCutPoints--;
6671 m_log <<
"removing two cut points, cell " << cellId <<
" edge " << edge << endl;
6672 m_log << a_coordinate(cellId, 0) <<
" " << a_coordinate(cellId, 1) <<
" ";
6673 IF_CONSTEXPR(nDim == 3)
m_log << a_coordinate(cellId, 2) << endl;
6675 candidates[cndt].noCutPointsOnEdge[edge]++;
6682 MBool newCutPoint =
true;
6683 for(
MInt i = 0; i < candidates[cndt].noCutPointsOnEdge[edge]; i++) {
6684 MBool equal = (
ABS(pP[0] - cutPointsEdge(i, 0)) <= eps &&
ABS(pP[1] - cutPointsEdge(i, 1)) <= eps);
6685 IF_CONSTEXPR(nDim == 3) equal = (equal && (
ABS(pP[2] - cutPointsEdge(i, 2)) <= eps));
6687 newCutPoint =
false;
6693 const MInt cutPointNo = candidates[cndt].noCutPoints;
6694 candidates[cndt].cutEdges[cutPointNo] = edge;
6695 candidates[cndt].cutBodyIds[cutPointNo] = 0;
6696 for(
MInt i = 0; i < nDim; i++) {
6697 candidates[cndt].cutPoints[cutPointNo][i] = pP[i];
6699 candidates[cndt].noCutPoints++;
6700 candidates[cndt].noCutPointsOnEdge[edge]++;
6704 cerr <<
"** Warning fvbndrycndxd: " << endl;
6705 cerr <<
"cell " << cellId << endl;
6706 cerr <<
" -> Boundary candidate " << cndt <<
" has more than " << m_noEdges <<
" cut points" << endl;
6707 cerr <<
" -> Additional cut points are neglected" << endl;
6708 cerr << cellId <<
" " << a_coordinate(cellId, 0) <<
" " << a_coordinate(cellId, 1) <<
" ";
6709 IF_CONSTEXPR(nDim == 3) cerr << a_coordinate(cellId, 2) <<
" ";
6710 cerr << grid().tree().level(cellId) << endl;
6711 mTerm(1, AT_,
"Boundary cell has more than m_noEdges cut points");
6718 m_scaledCutCell =
true;
6726template <MInt nDim_>
6728 const MInt* candidateIds) {
6735 const MInt debugTimeStep = -2;
6736 const MInt debugGlobalId = 216221;
6738 for(
MInt cellId = 0; cellId < grid().tree().size(); cellId++) {
6739 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
6740 const MInt cand = candidateIds[cellId];
6742 cerr <<
" Original node-values: ";
6743 for(
MInt node = 0; node < m_noCorners; node++) {
6744 cerr << candidates[cand].nodalValues[0][node] <<
" ";
6750 std::ignore = candidateIds[0];
6753 ASSERT(nDim == 3 && nDim_ == 3,
"Intended otherwise!");
6757 const constexpr MInt childParentNode[8] = {0, 1, 2, 3, 4, 5, 6, 7};
6759 const constexpr MInt faceCornerCode[6][4] = {{0, 2, 6, 4}, {3, 1, 5, 7}, {1, 0, 4, 5},
6760 {2, 3, 7, 6}, {0, 1, 3, 2}, {5, 4, 6, 7}};
6763 const constexpr MInt noCorrectNodes = nDim == 2 ? 2 : 4;
6766 const constexpr MInt edge2Corner[12][2] = {{0, 2}, {1, 3}, {0, 1}, {2, 3}, {4, 6}, {5, 7},
6767 {4, 5}, {6, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
6770 const constexpr MInt corner2Edge[3][8] = {
6771 {0, 1, 0, 1, 4, 5, 4, 5}, {2, 2, 3, 3, 6, 6, 7, 7}, {8, 9, 10, 11, 8, 9, 10, 11}};
6774 for(
MInt id = 0;
id < (
MInt)m_cutLvlJumpCandidates.size();
id++) {
6775 const MInt parentCand = m_cutLvlJumpCandidates[
id].parentCandId;
6776 const MInt candId = m_cutLvlJumpCandidates[
id].candId;
6777 const MInt cellId = candidates[candId].cellId;
6778 ASSERT(grid().tree().parent(cellId) == candidates[parentCand].cellId,
"");
6779 ASSERT(candidates[parentCand].isbndryLvlJumpParent,
"");
6781 for(
MInt node = 0; node < m_noCorners; node++) {
6782 ASSERT(candidates[parentCand].nodeValueSet[node],
"");
6785 const MInt noJumps = m_cutLvlJumpCandidates[
id].noJumps;
6786 ASSERT(noJumps <= 7 && noJumps > 0, to_string(noJumps));
6788 const MInt childPos = m_cutLvlJumpCandidates[
id].childId;
6791 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
6792 cerr <<
"Cell-Info: noJumps" << noJumps << endl;
6793 for(
MInt j = 0; j < noJumps; j++) {
6794 cerr <<
"dir " << m_cutLvlJumpCandidates[
id].dirs[j] <<
" type " << m_cutLvlJumpCandidates[
id].neighborType[j];
6801 const MInt parentNode = childParentNode[childPos];
6802 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6803 candidates[candId].nodalValues[set][parentNode] = candidates[parentCand].nodalValues[set][parentNode];
6806 const MInt parentEdgeIds[3] = {corner2Edge[0][parentNode], corner2Edge[1][parentNode], corner2Edge[2][parentNode]};
6809 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
6810 cerr <<
" Child-Pos " << childPos <<
" setting (1) node " << parentNode <<
" to "
6811 << candidates[parentCand].nodalValues[0][parentNode] << endl;
6816 for(
MInt j = 0; j < noJumps; j++) {
6817 const MInt jumpDir = m_cutLvlJumpCandidates[
id].dirs[j];
6818 const MInt type = m_cutLvlJumpCandidates[
id].neighborType[j];
6821 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
6822 cerr << childPos <<
" " << noJumps << type << endl;
6831 for(
MInt i = 0; i < noCorrectNodes; i++) {
6832 const MInt node = faceCornerCode[jumpDir][i];
6834 if(node == parentNode)
continue;
6836 MBool isDiagonalNode =
false;
6839 const MInt secondDir = m_cutLvlJumpCandidates[
id].diagonalDirs[j];
6840 for(
MInt l = 0; l < noCorrectNodes; l++) {
6841 if(node == faceCornerCode[secondDir][l]) {
6842 isDiagonalNode =
true;
6846 if(!isDiagonalNode)
continue;
6849 MInt edgeInterpolation = -1;
6850 for(
MInt k = 0; k < 3; k++) {
6851 MInt parentEdgeCorner = edge2Corner[parentEdgeIds[k]][0];
6852 if(parentEdgeCorner == parentNode) {
6853 parentEdgeCorner = edge2Corner[parentEdgeIds[k]][1];
6855 if(parentEdgeCorner == node) {
6856 edgeInterpolation = k;
6861 if(edgeInterpolation > -1) {
6862 const MInt parentEdgeId = parentEdgeIds[edgeInterpolation];
6863 const MInt parentNodes[2] = {edge2Corner[parentEdgeId][0], edge2Corner[parentEdgeId][1]};
6864 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6866 * (candidates[parentCand].nodalValues[set][parentNodes[0]]
6867 + candidates[parentCand].nodalValues[set][parentNodes[1]]);
6868 candidates[candId].nodalValues[set][node] = phi;
6872 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
6873 cerr <<
" (2 ) Setting node " << node <<
" to " << candidates[candId].nodalValues[0][node] << endl;
6878 ASSERT(type == 0,
"");
6879 const MInt parentNodes[4] = {faceCornerCode[jumpDir][0], faceCornerCode[jumpDir][1],
6880 faceCornerCode[jumpDir][2], faceCornerCode[jumpDir][3]};
6881 for(
MInt set = 0; set < m_noLevelSetsUsedForMb; set++) {
6883 * (candidates[parentCand].nodalValues[set][parentNodes[0]]
6884 + candidates[parentCand].nodalValues[set][parentNodes[1]]
6885 + candidates[parentCand].nodalValues[set][parentNodes[2]]
6886 + candidates[parentCand].nodalValues[set][parentNodes[3]]);
6887 candidates[candId].nodalValues[set][node] = phi;
6891 if(
globalTimeStep > debugTimeStep && grid().tree().globalId(cellId) == debugGlobalId) {
6892 cerr <<
" (3 ) Setting node " << node <<
" to " << candidates[candId].nodalValues[0][node] << endl;
6901 const MInt reverseNode[3][8] = {{1, 0, 3, 2, 5, 4, 7, 6}, {2, 3, 0, 1, 6, 7, 4, 5}, {4, 5, 6, 7, 0, 1, 2, 3}};
6905 std::set<std::pair<MInt, MInt>> nghbrSet;
6907 MInt nghbrIds[m_noCorners]{};
6908 MInt nghbrNodes[m_noCorners]{};
6910 const MInt nodeStencil[3][8] = {{0, 1, 0, 1, 0, 1, 0, 1}, {2, 2, 3, 3, 2, 2, 3, 3}, {4, 4, 4, 4, 5, 5, 5, 5}};
6913 for(
MInt cnd = 0; cnd < (signed)candidates.size(); cnd++) {
6914 const MInt cellId = candidates[cnd].cellId;
6915 for(
MInt node = 0; node < m_noCorners; node++) {
6917 nghbrSet.insert(std::make_pair(cellId, node));
6919 for(
MInt i = 0; i < nDim; i++) {
6920 const MInt firstDir = nodeStencil[i][node];
6921 const MInt firstNghbrNode = reverseNode[i][node];
6922 const MInt firstNghbrId = grid().tree().neighbor(cellId, firstDir);
6923 if(firstNghbrId > -1) {
6924 nghbrSet.insert(make_pair(firstNghbrId, firstNghbrNode));
6925 for(
MInt j = 0; j < nDim; j++) {
6926 const MInt secondDir = nodeStencil[j][node];
6927 if(secondDir == firstDir)
continue;
6928 const MInt secondNghbrNode = reverseNode[j][firstNghbrNode];
6929 const MInt secondNghbrId = grid().tree().neighbor(firstNghbrId, secondDir);
6930 if(secondNghbrId > -1) {
6931 nghbrSet.insert(make_pair(secondNghbrId, secondNghbrNode));
6932 IF_CONSTEXPR(nDim == 3) {
6933 for(
MInt k = 0; k < nDim; k++) {
6934 const MInt thirdDir = nodeStencil[k][node];
6935 if(thirdDir == firstDir || thirdDir == secondDir)
continue;
6936 const MInt thirdNghbrNode = reverseNode[k][secondNghbrNode];
6937 const MInt thirdNghbrId = grid().tree().neighbor(secondNghbrId, thirdDir);
6938 if(thirdNghbrId > -1) {
6939 nghbrSet.insert(make_pair(thirdNghbrId, thirdNghbrNode));
6951 MInt noNeighborsPerNode = 0;
6952 for(
auto it = nghbrSet.begin(); it != nghbrSet.end(); it++) {
6953 nghbrIds[noNeighborsPerNode] = it->first;
6954 nghbrNodes[noNeighborsPerNode] = it->second;
6955 noNeighborsPerNode++;
6958 for(
MInt nghbrNode = 0; nghbrNode < noNeighborsPerNode; nghbrNode++) {
6959 const MInt nghbrCand = candidateIds[nghbrIds[nghbrNode]];
6960 if(nghbrCand < 0)
continue;
6961 if(fabs(candidates[nghbrCand].nodalValues[0][nghbrNodes[nghbrNode]] - candidates[cnd].nodalValues[0][node])
6963 cerr <<
"Nodal-value missmatch " << cellId <<
" " << grid().tree().globalId(cellId) <<
" " << node <<
" "
6964 << candidates[cnd].nodalValues[0][node] <<
" " << nghbrIds[nghbrNode] <<
" "
6965 << grid().tree().globalId(nghbrIds[nghbrNode]) <<
" " << nghbrNodes[nghbrNode] <<
" "
6966 << candidates[nghbrCand].nodalValues[0][nghbrNodes[nghbrNode]] << endl;
6979template <MInt nDim_>
6981 MInt* nghbrIds,
MInt* nghbrNodes) {
6984 const MInt nodeStencil[3][8] = {{0, 1, 0, 1, 0, 1, 0, 1}, {2, 2, 3, 3, 2, 2, 3, 3}, {4, 4, 4, 4, 5, 5, 5, 5}};
6986 const MInt reverseNode[3][8] = {{1, 0, 3, 2, 5, 4, 7, 6}, {2, 3, 0, 1, 6, 7, 4, 5}, {4, 5, 6, 7, 0, 1, 2, 3}};
6988 std::set<std::pair<MInt, MInt>> nghbrSet;
6992 nghbrSet.insert(std::make_pair(cellId, node));
6994 for(
MInt i = 0; i < nDim; i++) {
6995 const MInt firstDir = nodeStencil[i][node];
6996 const MInt firstNghbrNode = reverseNode[i][node];
6997 const MInt firstNghbrId = grid().tree().neighbor(cellId, firstDir);
6998 if(firstNghbrId > -1) {
6999 nghbrSet.insert(make_pair(firstNghbrId, firstNghbrNode));
7000 for(
MInt j = 0; j < nDim; j++) {
7001 const MInt secondDir = nodeStencil[j][node];
7002 if(secondDir == firstDir)
continue;
7003 const MInt secondNghbrNode = reverseNode[j][firstNghbrNode];
7004 const MInt secondNghbrId = grid().tree().neighbor(firstNghbrId, secondDir);
7005 if(secondNghbrId > -1) {
7006 nghbrSet.insert(make_pair(secondNghbrId, secondNghbrNode));
7007 IF_CONSTEXPR(nDim == 3) {
7008 for(
MInt k = 0; k < nDim; k++) {
7009 const MInt thirdDir = nodeStencil[k][node];
7010 if(thirdDir == firstDir || thirdDir == secondDir)
continue;
7011 const MInt thirdNghbrNode = reverseNode[k][secondNghbrNode];
7012 const MInt thirdNghbrId = grid().tree().neighbor(secondNghbrId, thirdDir);
7013 if(thirdNghbrId > -1) {
7014 nghbrSet.insert(make_pair(thirdNghbrId, thirdNghbrNode));
7025 noNeighborsPerNode = 0;
7026 for(
auto it = nghbrSet.begin(); it != nghbrSet.end(); it++) {
7027 nghbrIds[noNeighborsPerNode] = it->first;
7028 nghbrNodes[noNeighborsPerNode] = it->second;
7029 noNeighborsPerNode++;
std::array< MInt, maxNoCutPoints > cutBodyIds
std::array< MFloat, maxNoCartesianSurfaces > cartFaceArea
std::array< std::array< MInt, maxNoFaceVertices >, maxNoTotalSurfaces > allFacesPointIds
std::array< MInt, noFaces > noFacesPerCartesianDir
std::array< MInt, maxNoTotalSurfaces > allFacesNoPoints
std::array< MInt, maxNoCutPoints > cutEdges
std::array< MInt, maxNoBoundarySurfaces > boundarySurfaceBodyId
std::array< MFloat, nDim > volumetricCentroid
std::array< MBool, noFaces > externalFaces
std::array< MInt, maxNoBoundarySurfaces > boundarySurfaceBndryCndId
std::array< std::array< MFloat, nDim >, maxNoBoundarySurfaces > boundarySurfaceCentroid
std::array< std::array< MFloat, nDim >, maxNoBoundarySurfaces > boundarySurfaceNormal
std::array< std::array< MFloat, nDim >, maxNoCartesianSurfaces > cartFaceCentroid
std::array< MInt, maxNoCartesianSurfaces > cartFaceDir
std::array< MFloat, maxNoBoundarySurfaces > boundarySurfaceArea
std::array< std::array< MBool, noCorners >, maxNoSets > cornerIsInsideGeometry
std::array< MInt, maxNoSets > associatedBodyIds
std::array< MInt, maxNoTotalSurfaces > allFacesBodyId
std::array< std::array< MFloat, nDim >, maxNoCutPoints > cutPoints
std::vector< MInt > vertices
std::array< MFloat, nDim > center
std::array< MFloat, nDim > normal
void compVolumeIntegrals_pyraBased3(std::vector< polyCutCell > *, std::vector< polyFace > *, const std::vector< polyVertex > *)
void writeVTKFileOfCell(MInt cellId, std::vector< polyFace > *faces, const std::vector< polyVertex > *vertices, MInt set)
void addPoint(const std::vector< polyVertex > *, const MInt *, const MInt, std::array< MFloat, nDim >)
adds an additional MC point in the cell sums up all n points passed in indices. points have to be sto...
void writeInfo(std::vector< CutCell< nDim > > &cutCellData, MUint, MInt)
void computeCutFaces(std::vector< CutCell< nDim > > &cutCellData, const MInt maxNoSurfaces, const MInt tCutGroup)
void exchangeNodalValues(const MInt **maxLevelWindowCells, const MInt *noMaxLevelWindowCells, const MInt **maxLevelHaloCells, std::vector< CutCandidate< nDim > > &candidates, MInt *candidateIds)
halo cell - window cell exchange of nodal scalar field data
void fillCutCellData(std::vector< CutCandidate< nDim > > &candidates, std::vector< CutCell< nDim > > &cutCellData, std::vector< MInt > cutCellIdMapping)
void getNeighborNodes(const MInt, const MInt, MInt, MInt *, MInt *)
gets all neighbor cellIds and matching nodes for a given cell and node
void computeNodalValues(std::vector< CutCandidate< nDim > > &candidates, MInt *candidateIds, const MFloat *scalarField, const MInt *bodyIdField, const MBool *const gapPropertyField, const MBool gapClosure)
1) add cutCandidates based on the solver m_bndryCandidates 2) set cutCandidate information 3) set the...
void compFaceIntegrals_pyraBased3(polyFace *face, const std::vector< polyVertex > *vertices, MFloat *MC, MFloat *VC, MFloat *XC)
void correctNodalValuesAtLevelJump(std::vector< CutCandidate< nDim > > &candidates, const MInt *)
corrects the nodal values at bndry level jumps for ls-based geometries
std::conditional< nDim_==3, polyCutCell3D, polyCutCell2D >::type polyCutCell
void computeCutPointsFromSTL(std::vector< CutCandidate< nDim > > &candidates)
computes cut points where candidate intersects with the geometry note: can not handle bndryRefinement...
void computeCutPoints(std::vector< CutCandidate< nDim > > &candidates, const MInt *candidateIds, std::vector< MInt > &candidatesOrder)
Computes cutpoints for the scalar-Field with grid edges and updates cutPoint-Info in the cutCandidate...
MBool computeCutFaceSimple(CutCell< nDim > &cutCell)
This class is a ScratchSpace.
void fill(T val)
fill the scratch with a given value
pointer p
Deprecated: use [] instead!
void mTerm(const MInt errorCode, const MString &location, const MString &message)
constexpr Real POW2(const Real x)
MBool approx(const T &, const U &, const T)
constexpr T mMin(const T &x, const T &y)
constexpr T mMax(const T &x, const T &y)
constexpr MLong IPOW2(MInt x)
int MPI_Isend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Isend
int MPI_Irecv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Irecv
int MPI_Waitall(int count, MPI_Request *request, MPI_Status *status, const MString &name)
same as MPI_Waitall
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
Namespace for auxiliary functions/classes.