computes the following boundary cell member variables: 2 cut points are assumed for each cell face multiple boundary cut faces are allowed
cell structure is determined by using a modified marching cubes approach (see e.g. Evgeni V. Chernyaev, Marching Cubes 33 and further publications) uses fvlookuptable.h
extended version of createCutFaceMG() method which is able to handle split cells and ambiguous cell configurations correctly.
24899 {
24900 TRACE();
24901
24902 static constexpr MInt cornerIndices[8][3] = {{-1, -1, -1}, {1, -1, -1}, {-1, 1, -1}, {1, 1, -1},
24903 {-1, -1, 1}, {1, -1, 1}, {-1, 1, 1}, {1, 1, 1}};
24904
24905 static constexpr MInt cornersMCtoSOLVER[8] = {2, 0, 1, 3, 6, 4, 5, 7};
24906 static constexpr MInt cornersSOLVERtoMC[8] = {1, 2, 0, 3, 5, 6, 4, 7};
24907 static constexpr MInt edgesMCtoSOLVER[12] = {0, 2, 1, 3, 4, 6, 5, 7, 10, 8, 9, 11};
24908 static constexpr MInt facesMCtoSOLVER[6] = {0, 2, 1, 3, 4, 5};
24909 static constexpr MInt neighborCorner[6] = {1, -1, 2, -2, 4, -4};
24910 static constexpr MInt cornerFaceMapping[24] = {0, 2, 4, 1, 2, 4, 0, 3, 4, 1, 3, 4,
24911 0, 2, 5, 1, 2, 5, 0, 3, 5, 1, 3, 5};
24912 static constexpr MInt faceCornerMapping[24] = {0, 2, 4, 6, 1, 3, 5, 7, 0, 1, 4, 5,
24913 2, 3, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7};
24914
24915 unsigned char outcode_MC = 0;
24919 MFloat corner[8][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
24921 MInt presentCases[15] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
24922 MInt* presentSubCases[15];
24923 const MInt maxNoSubCases = 48;
24924 MIntScratchSpace presentSubCases_scratch(15 * maxNoSubCases, AT_,
"presentSubCases_scratch");
24925 for(
MInt i = 0; i < 15; i++)
24926 presentSubCases[i] = &presentSubCases_scratch.p[maxNoSubCases * i];
24927 for(
MInt i = 0; i < 15; i++) {
24928 for(
MInt j = 0; j < maxNoSubCases; j++) {
24929 presentSubCases[i][j] = 0;
24930 }
24931 }
24932 MBool currentOutcode =
false;
24938 MFloat* faceCentroid[6];
24940 for(
MInt i = 0; i < 6; i++)
24941 faceCentroid[i] = &faceCentroid_scratch.p[i * nDim];
24942 MFloat* faceCentroid_0[6];
24944 for(
MInt i = 0; i < 6; i++)
24945 faceCentroid_0[i] = &faceCentroid_0_scratch.
p[i * nDim];
24946 MFloat* faceCentroid_00[6];
24947 MFloatScratchSpace faceCentroid_00_scratch(nDim * 6, AT_,
"faceCentroid_00_scratch");
24948 for(
MInt i = 0; i < 6; i++)
24949 faceCentroid_00[i] = &faceCentroid_00_scratch.
p[i * nDim];
24952 for(
MInt i = 0; i < 6; i++)
24953 normal[i] = &normal_scratch.
p[i * nDim];
24954 const MInt maxNoPyramids = 13;
24955 MFloat* pyraCentroid[maxNoPyramids];
24956 MFloatScratchSpace pyraCentroid_scratch(nDim * maxNoPyramids, AT_,
"pyraCentroid_scratch");
24957 for(
MInt i = 0; i < maxNoPyramids; i++)
24958 pyraCentroid[i] = &pyraCentroid_scratch.
p[i * nDim];
24962 MFloat* pyraVolume = pyraVolume_scratch.getPointer();
24964 MFloat* faceVolume_0 = faceVolume_0_scratch.getPointer();
24966 MFloat* faceVolume_00 = faceVolume_00_scratch.getPointer();
24969 MFloat* cellCentroid_0 = cellCentroid_0_scratch.getPointer();
24971 MInt currentSubCase;
24990
24991 MInt splitCorner1 = -1;
24992 MInt splitCorner2 = -1;
24993 MInt splitCorner12 = -1;
24994 MInt splitCorner22 = -1;
25002 MInt* facepointers[6] = {&face1, &face2, &face3, &face4, &face5, &face6};
25003 MInt cutPoints[12] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1};
25005 MFloat* normalVec_c = normalVec_c_scratch.getPointer();
25008 MFloat* coordinates_c = coordinates_c_scratch.getPointer();
25011 MFloat* coordinates_Cell = coordinatesCell_scratch.getPointer();
25012 MFloat M[3] = {0, 0, 0};
25013 MFloat dummy_1[3] = {0, 0, 0};
25014 MBool nfs_cur[6] = {
false,
false,
false,
false,
false,
false};
25015 MBool nfs_cur_0[6] = {
false,
false,
false,
false,
false,
false};
25016 MBool nfs_cur_00[6] = {
false,
false,
false,
false,
false,
false};
25017 MBool faceCut[6] = {
false,
false,
false,
false,
false,
false};
25018 MBool faceCut_0[6] = {
false,
false,
false,
false,
false,
false};
25019 MBool faceCut_00[6] = {
false,
false,
false,
false,
false,
false};
25023 MFloat* faceDiff = faceDiff_scratch.getPointer();
25024
25025
25026 static constexpr MInt opposite[6] = {1, 0, 3, 2, 5, 4};
25027 MInt bndryId2 = -1, bndryId3 = -1;
25028 MInt cellId2 = -1, cellId3 = -1;
25031 MIntScratchSpace disambiguation_scratch(noCells, AT_,
"disambiguation_scratch");
25033 MIntScratchSpace ambiguousCells_scratch(noCells, AT_,
"ambiguousCells_scratch");
25036 MInt* splitFace = splitFace_scratch.getPointer();
25037 MInt* splitCell = splitCell_scratch.getPointer();
25038 MInt* disambiguation = disambiguation_scratch.getPointer();
25039 MInt* ambIds = ambIds_scratch.getPointer();
25040 MInt* ambiguousCells = ambiguousCells_scratch.getPointer();
25041 MInt* subCases = subCases_scratch.getPointer();
25042 MInt* mcCases = mcCases_scratch.getPointer();
25044 ambiguousCells[i] = -1;
25045 disambiguation[i] = -1;
25046 ambIds[i] = -1;
25047 splitCell[i] = 0;
25048 splitFace[i] = 0;
25049 }
25050 MInt noSplitCells = 0;
25051 MInt noSplitFaces = 0;
25052 MInt noAmbiguousCells = 0;
25057 stack<MInt> ambStack;
25059 MFloat* splitFaceVolume = splitFaceVolume_scratch.getPointer();
25060 MFloat* splitFaceCentroid[2];
25061 MFloatScratchSpace splitFaceCentroid_scratch(nDim * 2, AT_,
"splitFaceCentroid_scratch");
25062 for(
MInt i = 0; i < 2; i++)
25063 splitFaceCentroid[i] = &splitFaceCentroid_scratch.p[i * nDim];
25064 MFloat* splitFaceNormal[2];
25065 MFloatScratchSpace splitFaceNormal_scratch(nDim * 2, AT_,
"splitFaceNormal_scratch");
25066 for(
MInt i = 0; i < 2; i++)
25067 splitFaceNormal[i] = &splitFaceNormal_scratch.
p[i * nDim];
25070 MFloat* splitFaceCentroid2[2];
25071 MFloatScratchSpace splitFaceCentroid2_scratch(nDim * 2, AT_,
"splitFaceCentroid2_scratch");
25072 for(
MInt i = 0; i < 2; i++)
25073 splitFaceCentroid2[i] = &splitFaceCentroid2_scratch.p[i * nDim];
25074 MFloat* splitFaceNormal2[2];
25075 MFloatScratchSpace splitFaceNormal2_scratch(nDim * 2, AT_,
"splitFaceNormal2_scratch");
25076 for(
MInt i = 0; i < 2; i++)
25077 splitFaceNormal2[i] = &splitFaceNormal2_scratch.
p[i * nDim];
25078 MInt splitFaceId = -1;
25079 MInt splitFaceId2 = -1;
25080 MInt splitSurfaceIndexCounter = 2;
25082 MBool keepNghbrs =
true;
25083 MInt disamb_tmp = 0;
25084 MInt disambPointer[16][21];
25085 const MInt* disambPointer_dummy =
nullptr;
25087
25088
25089
25090 for(
MInt bndryId = 0; bndryId < noCells; bndryId++) {
25093 outcode_MC = 0;
25094
25097 }
25098
25100 continue;
25101 }
25104 continue;
25105 }
25106
25107
25108
25109
25110
25111 for(
MInt i = 0; i < 8; i++) {
25112 for(
MInt dim = 0; dim < nDim; dim++) {
25113 corner[i][dim] =
25115 }
25117 if(currentOutcode) outcode_MC = outcode_MC | (1 << cornersSOLVERtoMC[i]);
25118 }
25119
25120
25121 currentCase = (
MInt)cases[outcode_MC][0];
25122 mcCases[bndryId] = currentCase;
25123 currentSubCase = (
MInt)cases[outcode_MC][1];
25124 subCases[bndryId] = currentSubCase;
25125 if(!caseStatesSTL[currentCase][0]) {
25126 cerr << "Error in createCutFaceMGC: Case not implemented: " << currentCase << endl;
25127 continue;
25128 }
25129
25130
25131 disambiguation[bndryId] = -1;
25132 if(!caseStatesSTL[mcCases[bndryId]][1]) {
25133 ambIds[bndryId] = noAmbiguousCells;
25134 ambiguousCells[noAmbiguousCells++] = bndryId;
25135 }
25136 }
25137
25138
25139
25140 MIntScratchSpace cornerCellMapping_scratch(8 * noAmbiguousCells, AT_,
"cornerCellMapping_scratch");
25141 MInt* cornerCellMapping = cornerCellMapping_scratch.getPointer();
25142 for(
MInt i = 0; i < noAmbiguousCells; i++) {
25143 bdAmbId = ambiguousCells[i];
25145 outcode_MC = 0;
25146
25147
25148 for(
MInt j = 0; j < 8; j++) {
25149 for(
MInt dim = 0; dim < nDim; dim++) {
25152 }
25154 if(currentOutcode) outcode_MC = outcode_MC | (1 << cornersSOLVERtoMC[j]);
25155 }
25156 for(
MInt j = 0; j < 8; j++) {
25157 mcCorner = cornersSOLVERtoMC[j];
25158 if(!(outcode_MC & 1 << mcCorner))
25159 cornerCellMapping[i * 8 + j] = -1;
25160 else
25161 cornerCellMapping[i * 8 + j] =
cellId;
25162 }
25163 }
25164
25165
25166 for(
MInt ambId = 0; ambId < noAmbiguousCells; ambId++) {
25167 bdAmbId = ambiguousCells[ambId];
25168 disambiguation[bdAmbId] = 0;
25170 currentCase = mcCases[bdAmbId];
25171 currentSubCase = subCases[bdAmbId];
25172
25174 switch(currentCase) {
25175 case 3:
25176 faceMax = 1;
25177 break;
25178 case 4:
25180 continue;
25181 case 6:
25182 faceMax = 1;
25183 break;
25184 case 7:
25185 faceMax = 3;
25186 break;
25187 case 10:
25188 faceMax = 2;
25189 break;
25190 case 12:
25191 faceMax = 2;
25192 break;
25193 default: {
25194 stringstream errorMessage;
25195 errorMessage << " Error in FvBndryCnd3D::createCutFaceMGC: cell in ambigous cells does not have an "
25196 "ambiguous case... cellId: "
25197 <<
cellId <<
" case: " << currentCase <<
" bndryId: " << bdAmbId << endl;
25199 mTerm(1, AT_, errorMessage.str());
25200 }
25201 }
25202
25203 for(
MInt i = 0; i < faceMax; i++) {
25204
25205
25206 switch(currentCase) {
25207 case 3:
25208 faceTMP = facesMCtoSOLVER[tiling3_1STL[currentSubCase][2 + i]];
25209 break;
25210
25211 case 6:
25212 faceTMP = facesMCtoSOLVER[tiling6_1STL[currentSubCase][2 + i]];
25213 break;
25214 case 7:
25215 faceTMP = facesMCtoSOLVER[tiling7_1STL[currentSubCase][3 + i]];
25216 break;
25217
25218 case 10:
25219 faceTMP = facesMCtoSOLVER[tiling10_1STL[currentSubCase][2 + i]];
25220 break;
25221
25222 case 12:
25223 faceTMP = facesMCtoSOLVER[tiling12_1STL[currentSubCase][2 + i]];
25224 break;
25225 default: {
25226 stringstream errorMessage;
25227 errorMessage << "ERROR: Switch variable 'currentCase' with value " << currentCase << " not matching any case."
25228 << endl;
25229 mTerm(1, AT_, errorMessage.str());
25230 }
25231 }
25233 if(nghbrCellId < 0 || !m_solver->a_hasNeighbor(cellId, faceTMP)) {
25234 stringstream errorMessage;
25235 errorMessage << " Error in FvBndryCnd3D::createCutFaceMGC: no neighbor in ambiguous direction... exiting."
25236 << endl;
25237 errorMessage <<
" cellId: " <<
cellId <<
", ambiguous face: " << faceTMP << endl;
25239 mTerm(1, AT_, errorMessage.str());
25240 }
25241
25242 for(
MInt dim = 0; dim < nDim; dim++) {
25244 }
25246
25247
25248 if(currentOutcode) {
25249 disambiguation[bdAmbId] +=
IPOW2(faceMax - i - 1);
25250 }
25251 }
25252 }
25253
25254
25255 for(
MInt ambId = 0; ambId < noAmbiguousCells; ambId++) {
25256 bdAmbId = ambiguousCells[ambId];
25258 currentCase = mcCases[bdAmbId];
25259 currentSubCase = subCases[bdAmbId];
25260
25261 switch(currentCase) {
25262 case 3:
25263 if(currentSubCase < 12 && disambiguation[bdAmbId] == 1) {
25264 splitCell[bdAmbId] = true;
25265 noSplitCells++;
25266 } else if(currentSubCase >= 12 && disambiguation[bdAmbId] == 1) {
25267 splitFace[bdAmbId] = true;
25268 noSplitFaces++;
25269 }
25270 break;
25271 case 4:
25272 if(currentSubCase < 4) {
25273 splitCell[bdAmbId] = true;
25274 noSplitCells++;
25275 }
25276 break;
25277 case 6:
25278 if(currentSubCase < 24 && disambiguation[bdAmbId] == 1) {
25279 splitCell[bdAmbId] = true;
25280 noSplitCells++;
25281 } else if(currentSubCase >= 24 && disambiguation[bdAmbId] == 1) {
25282 splitFace[bdAmbId] = true;
25283 noSplitFaces++;
25284 }
25285 break;
25286 case 7:
25287 if(currentSubCase < 8) {
25288 switch(disambiguation[bdAmbId]) {
25289 case 0:
25290
25291 break;
25292 case 1:
25293 case 2:
25294 case 4:
25295 splitFace[bdAmbId] = true;
25296 noSplitFaces++;
25297 break;
25298 case 3:
25299 case 5:
25300 case 6:
25301 splitCell[bdAmbId] = true;
25302 noSplitCells++;
25303 break;
25304 case 7:
25305 splitCell[bdAmbId] = true;
25306 noSplitCells++;
25307 noSplitCells++;
25308 break;
25309 default: {
25310 stringstream errorMessage;
25311 errorMessage << "ERROR: Switch variable 'disambiguation[bdAmbId]' with value " << disambiguation[bdAmbId]
25312 << " not matching any case." << endl;
25313 mTerm(1, AT_, errorMessage.str());
25314 }
25315 }
25316 } else {
25317 switch(disambiguation[bdAmbId]) {
25318 case 0:
25319
25320 break;
25321 case 1:
25322 case 2:
25323 case 4:
25324 splitFace[bdAmbId] = true;
25325 noSplitFaces++;
25326 break;
25327 case 3:
25328 case 5:
25329 case 6:
25330
25331 splitFace[bdAmbId] = true;
25332 noSplitFaces++;
25333 noSplitFaces++;
25334 break;
25335 case 7:
25336 splitCell[bdAmbId] = true;
25337 noSplitCells++;
25338 break;
25339 default: {
25340 stringstream errorMessage;
25341 errorMessage << "ERROR: Switch variable 'disambiguation[bdAmbId]' with value " << disambiguation[bdAmbId]
25342 << " not matching any case." << endl;
25343 mTerm(1, AT_, errorMessage.str());
25344 }
25345 }
25346 }
25347 break;
25348 case 10:
25349 if(disambiguation[bdAmbId] == 3) {
25350 splitCell[bdAmbId] = true;
25351 noSplitCells++;
25352 }
25353 if(disambiguation[bdAmbId] == 2 || disambiguation[bdAmbId] == 1) {
25354 stringstream errorMessage;
25355 errorMessage << " found case 10 cell with disambiguation " << disambiguation[bdAmbId] << " and cellId "
25356 <<
cellId <<
". This means cell is a case 10.2 cell which is not implemented. Exiting..."
25357 << endl;
25358 mTerm(1, AT_, errorMessage.str());
25359 }
25360 break;
25361 case 12:
25362 if(disambiguation[bdAmbId] == 3) {
25363 splitCell[bdAmbId] = true;
25364 noSplitCells++;
25365 }
25366 if(disambiguation[bdAmbId] == 2 || disambiguation[bdAmbId] == 1) {
25367 stringstream errorMessage;
25368 errorMessage << " found case 12 cell with disambiguation " << disambiguation[bdAmbId] << " and cellId "
25369 <<
cellId <<
". This means cell is a case 12.2 cell which is not implemented. Exiting..."
25370 << endl;
25371 mTerm(1, AT_, errorMessage.str());
25372 }
25373 break;
25374 default: {
25375 stringstream errorMessage;
25376 errorMessage << "ERROR: Switch variable 'currentCase' with value " << currentCase << " not matching any case."
25377 << endl;
25378 mTerm(1, AT_, errorMessage.str());
25379 }
25380 }
25381 }
25382
25383
25384
25385
25386
25387
25388
25389
25390
25391
25392
25393
25394 const MInt maxNoSplitChildrenPerCell = 3;
25395 mAlloc(
m_splitSurfaces, (2 + noSplitFaces) * 2 * maxNoSplitChildrenPerCell,
"m_splitSurfaces", -1, AT_);
25397 mAlloc(
m_splitChildren, (noCells + noSplitCells) * maxNoSplitChildrenPerCell,
"m_splitChildren", -1, AT_);
25398
25399
25400 for(
MInt bndryId = 0; bndryId < noCells; bndryId++) {
25401 noFaces = 0;
25404 outcode_MC = 0;
25405 volume_C = 0;
25406 coordinates_Cell[0] = 0;
25407 coordinates_Cell[1] = 0;
25408 coordinates_Cell[2] = 0;
25409
25410
25411
25412 for(
MInt cutPoint = 0; cutPoint <
m_bndryCells->
a[bndryId].m_srfcs[0]->m_noCutPoints; cutPoint++) {
25413 cutPoints[
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutEdge[cutPoint]] = cutPoint;
25414 }
25415
25416 for(
MInt i = 0; i < nDim; i++) {
25418 faceVolume[2 * i + 1] = faceVolume[2 * i];
25419 for(
MInt j = 0; j < 3; j++) {
25422 }
25425 }
25426
25427
25430 for(
MInt dim = 0; dim < 3; dim++) {
25432 }
25434 m_bndryCells->
a[bndryId].m_srfcs[0]->m_area = faceVolume[0];
25435 for(
MInt dim = 0; dim < 3; dim++) {
25436 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = faceCentroid[0][dim];
25437 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = F1 / sqrt(F3);
25438 }
25439 for(
MInt i = 0; i < nDim; i++) {
25441 }
25442 continue;
25443 }
25444
25445
25446 for(
MInt i = 0; i < 8; i++) {
25447 for(
MInt dim = 0; dim < nDim; dim++) {
25450 }
25451 }
25452 for(
MInt i = 0; i < 6; i++) {
25453 faceCut[i] = false;
25454 faceCut_0[i] = false;
25455 }
25456
25457
25458
25459 currentCase = mcCases[bndryId];
25460 presentCases[currentCase]++;
25461 currentSubCase = subCases[bndryId];
25462 presentSubCases[currentCase][currentSubCase]++;
25463
25464 if(!(
m_bndryCells->
a[bndryId].m_srfcs[0]->m_noCutPoints == caseCutPoints[currentCase])) {
25465 cerr << "FvBndryCnd3D::createCutFaceMGC - Error: number of cut points not correct!" << endl;
25466 cerr << "Case: " << currentCase << " cutPointsTheory: " << caseCutPoints[currentCase]
25467 <<
" actual Cut Points: " <<
m_bndryCells->
a[bndryId].m_srfcs[0]->m_noCutPoints << endl;
25469 stringstream fileName;
25470 fileName <<
cellId <<
".stl";
25472 cerr << "Wrote stl-file of cell. FileName: " << (fileName.str()) << endl;
25473 }
25474 if(!caseStatesSTL[currentCase][0]) {
25475 cerr << "FvBndryCnd3D::createCutFaceMGC - Error: Case not implemented: " << currentCase << endl;
25477 stringstream fileName;
25478 fileName <<
cellId <<
".stl";
25480 cerr << "Wrote stl-file of cell. FileName: " << (fileName.str()) << endl;
25481 continue;
25482 }
25483
25484
25485
25486 if(!caseStatesSTL[currentCase][1]) {
25487 switch(currentCase) {
25488 case 3: {
25489 for(
MInt i = 0; i < 6; i++) {
25490 faceVolume_0[i] = faceVolume[i];
25491 for(
MInt j = 0; j < 3; j++) {
25492 faceCentroid_0[i][j] = faceCentroid[i][j];
25493 }
25494 }
25495 cellVolume_0 = gridCellVolume;
25496 for(
MInt i = 0; i < 3; i++) {
25498 }
25499
25500
25501 if((currentSubCase < 12 && disambiguation[bndryId] == 0)
25502 || (currentSubCase >= 12 && disambiguation[bndryId] == 1)) {
25503
25504
25505
25507
25508 for(
MInt face = 0; face < 6; face++) {
25509 nfs_cur[face] = nfs3_2[currentSubCase][face];
25510 }
25511 noFaces = 5;
25512 p_0 = corner[cornersMCtoSOLVER[tiling3_2STL[currentSubCase][0]]];
25513 p_0s = corner[cornersMCtoSOLVER[tiling3_2STL[currentSubCase][1]]];
25514
25515 cutDummy = edgesMCtoSOLVER[tiling3_2STL[currentSubCase][2]];
25516 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25517 cutDummy = edgesMCtoSOLVER[tiling3_2STL[currentSubCase][3]];
25518 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25519 cutDummy = edgesMCtoSOLVER[tiling3_2STL[currentSubCase][4]];
25520 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25521 cutDummy = edgesMCtoSOLVER[tiling3_2STL[currentSubCase][5]];
25522 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25523 cutDummy = edgesMCtoSOLVER[tiling3_2STL[currentSubCase][6]];
25524 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25525 cutDummy = edgesMCtoSOLVER[tiling3_2STL[currentSubCase][7]];
25526 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25527 face1 = facesMCtoSOLVER[tiling3_2STL[currentSubCase][8]];
25528 face2 = facesMCtoSOLVER[tiling3_2STL[currentSubCase][9]];
25529 face3 = facesMCtoSOLVER[tiling3_2STL[currentSubCase][10]];
25530 face4 = facesMCtoSOLVER[tiling3_2STL[currentSubCase][11]];
25531 face5 = facesMCtoSOLVER[tiling3_2STL[currentSubCase][12]];
25532
25533 computeTri(p_0, p_1, p_2, &faceVolume[face1], faceCentroid[face1], normal[face1]);
25534 computeTri(p_0, p_2, p_3, &faceVolume[face2], faceCentroid[face2], normal[face2]);
25535 computeTri(p_0s, p_2s, p_3s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
25536 computeTri(p_0s, p_1s, p_2s, &faceVolume[face4], faceCentroid[face4], normal[face4]);
25537 computePoly6(p_0, p_1, p_3s, p_0s, p_1s, p_3, &faceVolume[face5], faceCentroid[face5], normal[face5]);
25538
25539 computePoly6(p_1, p_2, p_3, p_1s, p_2s, p_3s, &area_c, coordinates_c, normalVec_c);
25540
25541 maia::math::vecAvg<3>(M, p_0, p_0s, p_1, p_1s, p_2, p_2s, p_3, p_3s);
25542
25543 for(
MInt i = 0; i < 5; i++) {
25544 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
25545 &pyraVolume[i], pyraCentroid[i]);
25546 }
25547 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
25548
25549 for(
MInt i = 0; i < 6; i++) {
25550 volume_C += pyraVolume[i];
25551 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
25552 }
25553 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
25554
25555 if(currentSubCase >= 12) {
25556
25557 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
25558 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
25559 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
25560 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
25561 p_1ss = corner[cornersMCtoSOLVER[tiling3_2STL[currentSubCase][13]]];
25562 p_2ss = corner[cornersMCtoSOLVER[tiling3_2STL[currentSubCase][14]]];
25563 splitCorner1 = cornersMCtoSOLVER[tiling3_2STL[currentSubCase][13]];
25564 splitCorner2 = cornersMCtoSOLVER[tiling3_2STL[currentSubCase][14]];
25565 computeTri(p_1s, p_1ss, p_3, &splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0]);
25566 computeTri(p_1, p_2ss, p_3s, &splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1]);
25567 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
25568
25570 splitFaceId = face5;
25571 }
25572
25574
25576 for(
MInt dim = 0; dim < 3; dim++) {
25577 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
25578 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
25579 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
25580 }
25581
25582 faceCut[face1] = true;
25583 faceCut[face2] = true;
25584 faceCut[face3] = true;
25585 faceCut[face4] = true;
25586 faceCut[face5] = true;
25587
25588 absA = F0;
25589 for(
MInt i = 0; i < nDim; i++) {
25590 if(!nfs_cur[2 * i + 1]) faceVolume[2 * i + 1] = 0;
25591 if(!nfs_cur[2 * i]) faceVolume[2 * i] = 0;
25592 faceDiff[i] = faceVolume[2 * i + 1] - faceVolume[2 * i];
25593 absA +=
POW2(faceDiff[i]);
25594 }
25595
25596 if(currentSubCase < 12) {
25597 for(
MInt i = 0; i < nDim; i++) {
25598 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[i] = faceDiff[i] / sqrt(absA);
25599 }
25601 }
25602
25603 } else {
25604
25605 if(currentSubCase < 12) {
25606
25611 if(bndryId == 0)
25612 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! " << endl;
25615 } else {
25617 }
25618
25619 for(
MInt face = 0; face < 6; face++) {
25620 nfs_cur[face] = nfs3_1[currentSubCase][face];
25621 }
25622
25623
25624 subCaseDummy = tiling3_1STL[currentSubCase][0];
25625 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
25626 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
25627 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25628 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
25629 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25630 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
25631 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25632 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
25633 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
25634 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
25635
25636 faceCut[face1] = true;
25637 faceCut[face2] = true;
25638 faceCut[face3] = true;
25639
25640 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
25641 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
25642 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
25643
25644 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25645
25646 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25647
25648 if(currentSubCase >= 12) {
25649 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
25650 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
25651 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
25652
25653 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
25654
25656 } else {
25657 for(
MInt face = 0; face < 6; face++) {
25658 nfs_cur[face] = nfs1[subCaseDummy][face];
25659 }
25660 }
25661
25662
25664 for(
MInt dim = 0; dim < 3; dim++) {
25665 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
25666 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
25667 }
25669 for(
MInt dim = 0; dim < 3; dim++) {
25670 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
25671 }
25672
25673 if(currentSubCase >= 12) {
25674 faceVolume_0[face1] = faceVolume[face1];
25675 faceVolume_0[face2] = faceVolume[face2];
25676 faceVolume_0[face3] = faceVolume[face3];
25677 cellVolume_0 = volume_C;
25678 for(
MInt i = 0; i < 3; i++) {
25679 cellCentroid_0[i] = coordinates_Cell[i];
25680 faceCentroid_0[face1][i] = faceCentroid[face1][i];
25681 faceCentroid_0[face2][i] = faceCentroid[face2][i];
25682 faceCentroid_0[face3][i] = faceCentroid[face3][i];
25683 }
25684 }
25685
25686
25687 subCaseDummy = tiling3_1STL[currentSubCase][1];
25688 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
25689 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
25690 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25691 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
25692 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25693 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
25694 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25695 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
25696 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
25697 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
25698
25699 if(currentSubCase < 12) {
25700 for(
MInt face = 0; face < 6; face++) {
25701 nfs_cur_0[face] = nfs1[subCaseDummy][face];
25702 }
25703
25704 faceCut_0[face1] = true;
25705 faceCut_0[face2] = true;
25706 faceCut_0[face3] = true;
25707
25708 computeTri(p_0, p_1, p_3, &faceVolume_0[face1], faceCentroid_0[face1]);
25709 computeTri(p_0, p_3, p_2, &faceVolume_0[face2], faceCentroid_0[face2]);
25710 computeTri(p_0, p_2, p_1, &faceVolume_0[face3], faceCentroid_0[face3]);
25711
25712 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25713
25714 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25715
25716
25718 for(
MInt dim = 0; dim < 3; dim++) {
25719 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
25720 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
25721 }
25722
25724 for(
MInt dim = 0; dim < 3; dim++) {
25725 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
25726 }
25727
25728
25729 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]] = cellId2;
25730
25731 } else {
25732 faceCut[face1] = true;
25733 faceCut[face2] = true;
25734 faceCut[face3] = true;
25735
25736 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
25737 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
25738 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
25739
25740 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25741
25742 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25743
25744 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
25745 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
25746 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
25747
25748 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
25749
25751
25752
25755 for(
MInt dim = 0; dim < 3; dim++) {
25756 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
25757 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
25758 }
25759
25761 for(
MInt dim = 0; dim < 3; dim++) {
25762 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
25763 }
25764 }
25765 }
25766 break;
25767 }
25768 case 4: {
25769 for(
MInt i = 0; i < 6; i++) {
25770 faceVolume_0[i] = faceVolume[i];
25771 for(
MInt j = 0; j < 3; j++) {
25772 faceCentroid_0[i][j] = faceCentroid[i][j];
25773 }
25774 }
25775 cellVolume_0 = gridCellVolume;
25776 for(
MInt i = 0; i < 3; i++) {
25778 }
25779
25780
25781 if(currentSubCase < 4) {
25782
25787 if(bndryId == 0)
25788 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! " << endl;
25789
25790
25792
25793
25794 subCaseDummy = tiling4_1STL[currentSubCase][0];
25795 for(
MInt face = 0; face < 6; face++) {
25796 nfs_cur[face] = nfs1[subCaseDummy][face];
25797 }
25798 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
25799 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
25800 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25801 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
25802 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25803 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
25804 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25805 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
25806 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
25807 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
25808
25809 faceCut[face1] = true;
25810 faceCut[face2] = true;
25811 faceCut[face3] = true;
25812
25813 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
25814 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
25815 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
25816
25817 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25818
25819 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25820
25822 for(
MInt dim = 0; dim < 3; dim++) {
25823 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
25824 }
25825
25826
25828 for(
MInt dim = 0; dim < 3; dim++) {
25829 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
25830 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
25831 }
25832
25833
25834
25836
25837 subCaseDummy = tiling4_1STL[currentSubCase][1];
25838 for(
MInt face = 0; face < 6; face++) {
25839 nfs_cur_0[face] = nfs1[subCaseDummy][face];
25840 }
25841 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
25842 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
25843 p_1 =
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25844 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
25845 p_2 =
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25846 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
25847 p_3 =
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25848 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
25849 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
25850 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
25851
25852 faceCut_0[face1] = true;
25853 faceCut_0[face2] = true;
25854 faceCut_0[face3] = true;
25855
25856 computeTri(p_0, p_1, p_3, &faceVolume_0[face1], faceCentroid_0[face1]);
25857 computeTri(p_0, p_3, p_2, &faceVolume_0[face2], faceCentroid_0[face2]);
25858 computeTri(p_0, p_2, p_1, &faceVolume_0[face3], faceCentroid_0[face3]);
25859
25860 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25861
25862 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25863
25865 for(
MInt dim = 0; dim < 3; dim++) {
25866 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
25867 }
25868
25869
25871 for(
MInt dim = 0; dim < 3; dim++) {
25872 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
25873 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
25874 }
25875
25876
25877 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]] = cellId2;
25878
25879 } else {
25881 for(
MInt face = 0; face < 6; face++) {
25882 nfs_cur[face] = nfs4_1[currentSubCase][face];
25883 }
25884
25885 subCaseDummy = tiling4_1STL[currentSubCase][0];
25886 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
25887 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
25888 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25889 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
25890 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25891 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
25892 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25893 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
25894 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
25895 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
25896
25897 faceCut[face1] = true;
25898 faceCut[face2] = true;
25899 faceCut[face3] = true;
25900
25901 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
25902 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
25903 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
25904
25905 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25906
25907 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25908
25909 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
25910 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
25911 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
25912
25913 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
25914
25916
25917
25919 for(
MInt dim = 0; dim < 3; dim++) {
25920 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
25921 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
25922 }
25923
25924 faceVolume_0[face1] = faceVolume[face1];
25925 faceVolume_0[face2] = faceVolume[face2];
25926 faceVolume_0[face3] = faceVolume[face3];
25927 cellVolume_0 = volume_C;
25928 for(
MInt i = 0; i < 3; i++) {
25929 cellCentroid_0[i] = coordinates_Cell[i];
25930 faceCentroid_0[face1][i] = faceCentroid[face1][i];
25931 faceCentroid_0[face2][i] = faceCentroid[face2][i];
25932 faceCentroid_0[face3][i] = faceCentroid[face3][i];
25933 }
25934
25935
25936 subCaseDummy = tiling4_1STL[currentSubCase][1];
25937 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
25938 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
25939 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25940 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
25941 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25942 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
25943 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
25944 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
25945 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
25946 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
25947
25948 faceCut[face1] = true;
25949 faceCut[face2] = true;
25950 faceCut[face3] = true;
25951
25952 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
25953 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
25954 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
25955
25956 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
25957
25958 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
25959
25960 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
25961 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
25962 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
25963
25964 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
25965
25967
25968
25971 for(
MInt dim = 0; dim < 3; dim++) {
25972 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
25973 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
25974 }
25975
25977 for(
MInt dim = 0; dim < 3; dim++) {
25978 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
25979 }
25980 }
25981 break;
25982 }
25983 case 6: {
25984 for(
MInt i = 0; i < 6; i++) {
25985 faceVolume_0[i] = faceVolume[i];
25986 for(
MInt j = 0; j < 3; j++) {
25987 faceCentroid_0[i][j] = faceCentroid[i][j];
25988 }
25989 }
25990 cellVolume_0 = gridCellVolume;
25991 for(
MInt i = 0; i < 3; i++) {
25993 }
25994 if((currentSubCase < 24 && disambiguation[bndryId] == 1)
25995 || (currentSubCase >= 24 && disambiguation[bndryId] == 0)) {
25996
25997
25998 if(currentSubCase < 24) {
25999
26004 if(bndryId == 0)
26005 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! " << endl;
26008 } else {
26010 }
26011
26012 for(
MInt face = 0; face < 6; face++) {
26013 nfs_cur[face] = nfs6_1[currentSubCase][face];
26014 }
26015
26016
26017 subCaseDummy = tiling6_1STL[currentSubCase][0];
26018 p_0 = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][0]]];
26019 p_0s = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][1]]];
26020 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][2]];
26021 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26022 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][3]];
26023 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26024 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][4]];
26025 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26026 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][5]];
26027 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26028 face1 = facesMCtoSOLVER[tiling2STL[subCaseDummy][6]];
26029 face2 = facesMCtoSOLVER[tiling2STL[subCaseDummy][7]];
26030 face3 = facesMCtoSOLVER[tiling2STL[subCaseDummy][8]];
26031 face4 = facesMCtoSOLVER[tiling2STL[subCaseDummy][9]];
26032
26033 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2], normal[face2]);
26034 computeTri(p_0s, p_3s, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
26035 computeTrapez(p_0, p_0s, p_3s, p_3, &faceVolume[face1], faceCentroid[face1], normal[face1]);
26036 computeTrapez(p_2, p_2s, p_0s, p_0, &faceVolume[face4], faceCentroid[face4], normal[face4]);
26037
26038 computePoly4(p_3, p_3s, p_2s, p_2, &area_c, coordinates_c, normalVec_c);
26039
26040 maia::math::vecAvg<3>(M, p_0, p_0s, p_2, p_2s, p_3, p_3s);
26041
26042 for(
MInt i = 0; i < 4; i++) {
26043 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
26044 &pyraVolume[i], pyraCentroid[i]);
26045 }
26046 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[4], pyraCentroid[4]);
26047 for(
MInt i = 0; i < 5; i++) {
26048 volume_C += pyraVolume[i];
26049 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26050 }
26051 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26052
26053 if(currentSubCase >= 24) {
26054 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26055 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26056 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26057 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
26058
26059 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26060
26062 } else {
26063 for(
MInt face = 0; face < 6; face++) {
26064 nfs_cur[face] = nfs2[subCaseDummy][face];
26065 }
26066 }
26067
26068 faceCut[face1] = true;
26069 faceCut[face2] = true;
26070 faceCut[face3] = true;
26071 faceCut[face4] = true;
26072
26073
26075 for(
MInt dim = 0; dim < 3; dim++) {
26076 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26077 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26078 }
26080 for(
MInt dim = 0; dim < 3; dim++) {
26081 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26082 }
26083
26084 if(currentSubCase >= 24) {
26085 faceVolume_0[face1] = faceVolume[face1];
26086 faceVolume_0[face2] = faceVolume[face2];
26087 faceVolume_0[face3] = faceVolume[face3];
26088 faceVolume_0[face4] = faceVolume[face4];
26089 cellVolume_0 = volume_C;
26090 volume_C = 0;
26091 for(
MInt i = 0; i < 3; i++) {
26092 cellCentroid_0[i] = coordinates_Cell[i];
26093 faceCentroid_0[face1][i] = faceCentroid[face1][i];
26094 faceCentroid_0[face2][i] = faceCentroid[face2][i];
26095 faceCentroid_0[face3][i] = faceCentroid[face3][i];
26096 faceCentroid_0[face4][i] = faceCentroid[face4][i];
26097 coordinates_Cell[i] = 0;
26098 }
26099 }
26100
26101
26102 subCaseDummy = tiling6_1STL[currentSubCase][1];
26103 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
26104 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
26105 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26106 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
26107 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26108 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
26109 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26110 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
26111 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
26112 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
26113
26114 if(currentSubCase < 24) {
26115 for(
MInt face = 0; face < 6; face++) {
26116 nfs_cur_0[face] = nfs1[subCaseDummy][face];
26117 }
26118
26119 faceCut_0[face1] = true;
26120 faceCut_0[face2] = true;
26121 faceCut_0[face3] = true;
26122
26123 computeTri(p_0, p_1, p_3, &faceVolume_0[face1], faceCentroid_0[face1]);
26124 computeTri(p_0, p_3, p_2, &faceVolume_0[face2], faceCentroid_0[face2]);
26125 computeTri(p_0, p_2, p_1, &faceVolume_0[face3], faceCentroid_0[face3]);
26126
26127 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
26128
26129 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
26130
26131
26133 for(
MInt dim = 0; dim < 3; dim++) {
26134 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26135 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26136 }
26137
26139 for(
MInt dim = 0; dim < 3; dim++) {
26140 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
26141 }
26142
26143
26144 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]] = cellId2;
26145 } else {
26146 faceCut[face1] = true;
26147 faceCut[face2] = true;
26148 faceCut[face3] = true;
26149
26150 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
26151 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
26152 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
26153
26154 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
26155
26156 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
26157
26158 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26159 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26160 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26161
26162 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26163
26165
26166
26169 for(
MInt dim = 0; dim < 3; dim++) {
26170 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
26171 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
26172 }
26173
26175 for(
MInt dim = 0; dim < 3; dim++) {
26176 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26177 }
26178 }
26179
26180 } else {
26181
26182
26183
26184
26186
26187 for(
MInt face = 0; face < 6; face++) {
26188 nfs_cur[face] = nfs6_1[currentSubCase][face];
26189 }
26190 p_0 = corner[cornersMCtoSOLVER[tiling6_2STL[currentSubCase][0]]];
26191 p_1 = corner[cornersMCtoSOLVER[tiling6_2STL[currentSubCase][1]]];
26192 p_2 = corner[cornersMCtoSOLVER[tiling6_2STL[currentSubCase][2]]];
26193 p_4 = corner[cornersMCtoSOLVER[tiling6_2STL[currentSubCase][3]]];
26194 p_5 = corner[cornersMCtoSOLVER[tiling6_2STL[currentSubCase][4]]];
26195
26196 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][5]];
26197 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26198 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][6]];
26199 p_0ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26200 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][7]];
26201 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26202 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][8]];
26203 p_1ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26204 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][9]];
26205 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26206 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][10]];
26207 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26208 cutDummy = edgesMCtoSOLVER[tiling6_2STL[currentSubCase][11]];
26209 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26210 face1 = facesMCtoSOLVER[tiling6_2STL[currentSubCase][12]];
26211 face2 = facesMCtoSOLVER[tiling6_2STL[currentSubCase][13]];
26212 face3 = facesMCtoSOLVER[tiling6_2STL[currentSubCase][14]];
26213 face4 = facesMCtoSOLVER[tiling6_2STL[currentSubCase][15]];
26214 face5 = facesMCtoSOLVER[tiling6_2STL[currentSubCase][16]];
26215 face6 = facesMCtoSOLVER[tiling6_2STL[currentSubCase][17]];
26216
26217 computePoly6(p_1, p_1s, p_2s, p_2, p_2ss, p_1ss, &faceVolume[face1], faceCentroid[face1], normal[face1]);
26218 computeTrapez(p_0, p_1, p_1s, p_0s, &faceVolume[face2], faceCentroid[face2], normal[face2]);
26219 computeTrapez(p_0, p_1, p_1ss, p_0ss, &faceVolume[face3], faceCentroid[face3], normal[face3]);
26220 computeTri(p_2, p_3, p_2ss, &faceVolume[face4], faceCentroid[face4], normal[face4]);
26221 computeTri(p_2s, p_2, p_3, &faceVolume[face5], faceCentroid[face5], normal[face5]);
26222 computeTri(p_0, p_0s, p_0ss, &faceVolume[face6], faceCentroid[face6], normal[face6]);
26223
26224 maia::math::vecAvg<3>(M, p_0s, p_0ss, p_1s, p_1ss, p_2s, p_2ss, p_3, p_0, p_1, p_2);
26225
26226 for(
MInt i = 0; i < 6; i++) {
26227 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
26228 &pyraVolume[i], pyraCentroid[i]);
26229 }
26230
26231 computePoly4(p_0s, p_3, p_2s, p_1s, &area_c, coordinates_c, normalVec_c);
26232 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[6], pyraCentroid[6]);
26235 for(
MInt dim = 0; dim < 3; dim++) {
26236 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26237 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26238 }
26239
26240 computePoly4(p_0ss, p_1ss, p_2ss, p_3, &area_c, coordinates_c, normalVec_c);
26241 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[7], pyraCentroid[7]);
26245 for(
MInt dim = 0; dim < 3; dim++) {
26246 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
26247 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
26248 }
26249
26250 computePoly3(p_0s, p_0ss, p_3, &area_c, coordinates_c, normalVec_c);
26251 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[8], pyraCentroid[8]);
26255 for(
MInt dim = 0; dim < 3; dim++) {
26256 m_bndryCells->
a[bndryId].m_srfcs[2]->m_coordinates[dim] = coordinates_c[dim];
26257 m_bndryCells->
a[bndryId].m_srfcs[2]->m_normalVector[dim] = normalVec_c[dim];
26258 }
26259
26260 for(
MInt i = 0; i < 9; i++) {
26261 volume_C += pyraVolume[i];
26262 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26263 }
26264 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26265
26266 if(currentSubCase >= 24) {
26267
26268 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26269 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26270 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
26271 correctFace(&faceVolume[face5], faceCentroid[face5], &faceVolume_0[face5], faceCentroid_0[face5]);
26272 correctFace(&faceVolume[face6], faceCentroid[face6], &faceVolume_0[face6], faceCentroid_0[face6]);
26273 splitCorner1 = cornersMCtoSOLVER[tiling6_2STL[currentSubCase][3]];
26274 splitCorner2 = cornersMCtoSOLVER[tiling6_2STL[currentSubCase][4]];
26275 computeTri(p_4, p_2s, p_1s, &splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0]);
26276 computeTri(p_5, p_1ss, p_2ss, &splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1]);
26277
26278 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26279
26280 splitFaceId = face1;
26281 }
26282
26284 for(
MInt dim = 0; dim < 3; dim++) {
26285 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26286 }
26287
26288 faceCut[face1] = true;
26289 faceCut[face2] = true;
26290 faceCut[face3] = true;
26291 faceCut[face4] = true;
26292 faceCut[face5] = true;
26293 faceCut[face6] = true;
26294 }
26295
26296 break;
26297 }
26298 case 7: {
26299 disamb_tmp = disambiguation[bndryId];
26300 if(currentSubCase >= 8) disamb_tmp = 7 - disamb_tmp;
26301
26302 switch(disamb_tmp) {
26303 case 0: {
26304
26305
26306 for(
MInt i = 0; i < 6; i++) {
26307 faceVolume_0[i] = faceVolume[i];
26308 for(
MInt j = 0; j < 3; j++) {
26309 faceCentroid_0[i][j] = faceCentroid[i][j];
26310 }
26311 }
26312 cellVolume_0 = gridCellVolume;
26313 for(
MInt i = 0; i < 3; i++) {
26315 }
26316
26317 if(currentSubCase >= 8) {
26318
26323 if(bndryId == 0)
26324 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! "
26325 << endl;
26328 } else {
26330 }
26331
26332 for(
MInt face = 0; face < 6; face++) {
26333 nfs_cur[face] = nfs7_1[currentSubCase][face];
26334 }
26335
26336
26337 subCaseDummy = tiling7_4_1STL[currentSubCase][0];
26338 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
26339 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
26340 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26341 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
26342 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26343 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
26344 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26345 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
26346 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
26347 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
26348
26349 faceCut[face1] = true;
26350 faceCut[face2] = true;
26351 faceCut[face3] = true;
26352
26353 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
26354 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
26355 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
26356
26357 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
26358
26359 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
26360
26361 if(currentSubCase < 8) {
26362 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26363 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26364 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26365
26366 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26367
26369 } else {
26370 for(
MInt face = 0; face < 6; face++) {
26371 nfs_cur[face] = nfs1[subCaseDummy][face];
26372 }
26373 }
26374
26375
26377 for(
MInt dim = 0; dim < 3; dim++) {
26378 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26379 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26380 }
26382 for(
MInt dim = 0; dim < 3; dim++) {
26383 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26384 }
26385
26386 if(currentSubCase < 8) {
26387 faceVolume_0[face1] = faceVolume[face1];
26388 faceVolume_0[face2] = faceVolume[face2];
26389 faceVolume_0[face3] = faceVolume[face3];
26390 cellVolume_0 = volume_C;
26391 volume_C = 0;
26392 for(
MInt i = 0; i < 3; i++) {
26393 cellCentroid_0[i] = coordinates_Cell[i];
26394 faceCentroid_0[face1][i] = faceCentroid[face1][i];
26395 faceCentroid_0[face2][i] = faceCentroid[face2][i];
26396 faceCentroid_0[face3][i] = faceCentroid[face3][i];
26397 coordinates_Cell[i] = 0;
26398 }
26399 } else {
26400 volume_C = 0;
26401 for(
MInt i = 0; i < 3; i++) {
26402 coordinates_Cell[i] = 0;
26403 }
26404 }
26405
26406
26407 subCaseDummy = tiling7_4_1STL[currentSubCase][1];
26408 noFaces = 6;
26409 p_0 = corner[cornersMCtoSOLVER[tiling9STL[subCaseDummy][0]]];
26410 p_1 = corner[cornersMCtoSOLVER[tiling9STL[subCaseDummy][1]]];
26411 p_2 = corner[cornersMCtoSOLVER[tiling9STL[subCaseDummy][2]]];
26412 p_3 = corner[cornersMCtoSOLVER[tiling9STL[subCaseDummy][3]]];
26413
26414 cutDummy = edgesMCtoSOLVER[tiling9STL[subCaseDummy][4]];
26415 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26416 cutDummy = edgesMCtoSOLVER[tiling9STL[subCaseDummy][5]];
26417 p_1ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26418 cutDummy = edgesMCtoSOLVER[tiling9STL[subCaseDummy][6]];
26419 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26420 cutDummy = edgesMCtoSOLVER[tiling9STL[subCaseDummy][7]];
26421 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26422 cutDummy = edgesMCtoSOLVER[tiling9STL[subCaseDummy][8]];
26423 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26424 cutDummy = edgesMCtoSOLVER[tiling9STL[subCaseDummy][9]];
26425 p_3ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26426 face1 = facesMCtoSOLVER[tiling9STL[subCaseDummy][10]];
26427 face2 = facesMCtoSOLVER[tiling9STL[subCaseDummy][11]];
26428 face3 = facesMCtoSOLVER[tiling9STL[subCaseDummy][12]];
26429 face4 = facesMCtoSOLVER[tiling9STL[subCaseDummy][13]];
26430 face5 = facesMCtoSOLVER[tiling9STL[subCaseDummy][14]];
26431 face6 = facesMCtoSOLVER[tiling9STL[subCaseDummy][15]];
26432
26433 if(currentSubCase >= 8) {
26434 for(
MInt face = 0; face < 6; face++) {
26435 nfs_cur_0[face] = nfs9[subCaseDummy][face];
26436 faceCut_0[face] = true;
26437 }
26438
26439 computePoly5(p_0, p_3, p_3ss, p_2s, p_2, &faceVolume_0[face1], faceCentroid_0[face1], normal[face1]);
26440 computePoly5(p_0, p_1, p_1ss, p_3s, p_3, &faceVolume_0[face3], faceCentroid_0[face3], normal[face3]);
26441 computePoly5(p_0, p_2, p_2ss, p_1s, p_1, &faceVolume_0[face5], faceCentroid_0[face5], normal[face5]);
26442 computeTri(p_1, p_1s, p_1ss, &faceVolume_0[face2], faceCentroid_0[face2], normal[face2]);
26443 computeTri(p_2, p_2s, p_2ss, &faceVolume_0[face4], faceCentroid_0[face4], normal[face4]);
26444 computeTri(p_3, p_3s, p_3ss, &faceVolume_0[face6], faceCentroid_0[face6], normal[face6]);
26445
26446 computePoly6(p_1ss, p_1s, p_2ss, p_2s, p_3ss, p_3s, &area_c, coordinates_c, normalVec_c);
26447
26448 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_3, p_1s, p_2s, p_3s, p_1ss, p_2ss, p_3ss);
26449
26450 for(
MInt i = 0; i < 6; i++) {
26451 computePyra(&faceVolume_0[*facepointers[i]], faceCentroid_0[*facepointers[i]],
26452 normal[*facepointers[i]], M, &pyraVolume[i], pyraCentroid[i]);
26453 }
26454 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[6], pyraCentroid[6]);
26455 for(
MInt i = 0; i < 7; i++) {
26456 volume_C += pyraVolume[i];
26457 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26458 }
26459 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26460
26461
26463 for(
MInt dim = 0; dim < 3; dim++) {
26464 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26465 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26466 }
26467
26469 for(
MInt dim = 0; dim < 3; dim++) {
26470 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
26471 }
26472
26473
26474 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling9STL[subCaseDummy][0]]] = cellId2;
26475 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling9STL[subCaseDummy][1]]] = cellId2;
26476 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling9STL[subCaseDummy][2]]] = cellId2;
26477 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling9STL[subCaseDummy][3]]] = cellId2;
26478
26479 } else {
26480 for(
MInt face = 0; face < 6; face++) {
26481 faceCut[face] = true;
26482 }
26483
26484 computePoly5(p_0, p_3, p_3ss, p_2s, p_2, &faceVolume[face1], faceCentroid[face1], normal[face1]);
26485 computePoly5(p_0, p_1, p_1ss, p_3s, p_3, &faceVolume[face3], faceCentroid[face3], normal[face3]);
26486 computePoly5(p_0, p_2, p_2ss, p_1s, p_1, &faceVolume[face5], faceCentroid[face5], normal[face5]);
26487 computeTri(p_1, p_1s, p_1ss, &faceVolume[face2], faceCentroid[face2], normal[face2]);
26488 computeTri(p_2, p_2s, p_2ss, &faceVolume[face4], faceCentroid[face4], normal[face4]);
26489 computeTri(p_3, p_3s, p_3ss, &faceVolume[face6], faceCentroid[face6], normal[face6]);
26490
26491 computePoly6(p_1ss, p_1s, p_2ss, p_2s, p_3ss, p_3s, &area_c, coordinates_c, normalVec_c);
26492
26493 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_3, p_1s, p_2s, p_3s, p_1ss, p_2ss, p_3ss);
26494
26495 for(
MInt i = 0; i < 6; i++) {
26496 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]],
26497 M, &pyraVolume[i], pyraCentroid[i]);
26498 }
26499 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[6], pyraCentroid[6]);
26500 for(
MInt i = 0; i < 7; i++) {
26501 volume_C += pyraVolume[i];
26502 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26503 }
26504 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26505
26506 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26507 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26508 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26509 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
26510 correctFace(&faceVolume[face5], faceCentroid[face5], &faceVolume_0[face5], faceCentroid_0[face5]);
26511 correctFace(&faceVolume[face6], faceCentroid[face6], &faceVolume_0[face6], faceCentroid_0[face6]);
26512
26513 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26514
26516
26517
26520 for(
MInt dim = 0; dim < 3; dim++) {
26521 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
26522 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
26523 }
26524
26526 for(
MInt dim = 0; dim < 3; dim++) {
26527 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26528 }
26529 }
26530
26531 break;
26532 }
26533 case 1:
26534 case 2:
26535 case 4: {
26536 for(
MInt i = 0; i < 6; i++) {
26537 faceVolume_0[i] = faceVolume[i];
26538 faceVolume_00[i] = faceVolume[i];
26539 for(
MInt j = 0; j < 3; j++) {
26540 faceCentroid_0[i][j] = faceCentroid[i][j];
26541 }
26542 }
26543 cellVolume_0 = gridCellVolume;
26544 for(
MInt i = 0; i < 3; i++) {
26546 }
26547
26548 for(
MInt face = 0; face < 6; face++) {
26549 nfs_cur[face] = nfs7_1[currentSubCase][face];
26550 }
26551
26552 if(disamb_tmp == 1) disambPointer_dummy = &tiling7_3_1STL[0][0];
26553 if(disamb_tmp == 2) disambPointer_dummy = &tiling7_3_2STL[0][0];
26554 if(disamb_tmp == 4) disambPointer_dummy = &tiling7_3_4STL[0][0];
26555
26556 for(
MInt i = 0; i < 16; i++)
26557 for(
MInt j = 0; j < 21; j++)
26558 disambPointer[i][j] = disambPointer_dummy[i * 21 + j];
26559
26560 p_0 = corner[cornersMCtoSOLVER[disambPointer[currentSubCase][0]]];
26561 p_1 = corner[cornersMCtoSOLVER[disambPointer[currentSubCase][1]]];
26562 p_2 = corner[cornersMCtoSOLVER[disambPointer[currentSubCase][2]]];
26563 p_3 = corner[cornersMCtoSOLVER[disambPointer[currentSubCase][3]]];
26564 p_4 = corner[cornersMCtoSOLVER[disambPointer[currentSubCase][4]]];
26565 p_5 = corner[cornersMCtoSOLVER[disambPointer[currentSubCase][5]]];
26566 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][6]];
26567 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26568 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][7]];
26569 p_1ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26570 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][8]];
26571 p_1sss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26572 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][9]];
26573 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26574 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][10]];
26575 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26576 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][11]];
26577 p_2sss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26578 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][12]];
26579 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26580 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][13]];
26581 p_3ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26582 cutDummy = edgesMCtoSOLVER[disambPointer[currentSubCase][14]];
26583 p_3sss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26584
26585 face1 = facesMCtoSOLVER[disambPointer[currentSubCase][15]];
26586 face2 = facesMCtoSOLVER[disambPointer[currentSubCase][16]];
26587 face3 = facesMCtoSOLVER[disambPointer[currentSubCase][17]];
26588 face4 = facesMCtoSOLVER[disambPointer[currentSubCase][18]];
26589 face5 = facesMCtoSOLVER[disambPointer[currentSubCase][19]];
26590 face6 = facesMCtoSOLVER[disambPointer[currentSubCase][20]];
26591 if(currentSubCase < 8) {
26592
26594
26595
26596 computePoly6(p_1, p_1sss, p_3sss, p_3, p_3s, p_1ss, &faceVolume[face2], faceCentroid[face2],
26597 normal[face2]);
26598 computePoly6(p_1, p_1s, p_2ss, p_2, p_2sss, p_1sss, &faceVolume[face3], faceCentroid[face3],
26599 normal[face3]);
26600 computeTri(p_1, p_1ss, p_1s, &faceVolume[face4], faceCentroid[face4], normal[face4]);
26601 computeTri(p_2, p_2ss, p_2s, &faceVolume[face5], faceCentroid[face5], normal[face5]);
26602 computeTri(p_3, p_3ss, p_3s, &faceVolume[face6], faceCentroid[face6], normal[face6]);
26603 splitCorner1 = cornersMCtoSOLVER[disambPointer[currentSubCase][3]];
26604 splitCorner2 = cornersMCtoSOLVER[disambPointer[currentSubCase][2]];
26605 computeTri(p_3, p_3sss, p_3ss, &splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0]);
26606 computeTri(p_2, p_2s, p_2sss, &splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1]);
26607
26608 splitFaceId = face1;
26609 maia::math::vecAvg<3>(M, p_1s, p_1ss, p_1sss, p_2s, p_2ss, p_2sss, p_3s, p_3ss, p_3sss);
26610
26611 faceVolume[face1] = F0;
26612
26613 for(
MInt i = 0; i < 6; i++) {
26614 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]],
26615 M, &pyraVolume[i], pyraCentroid[i]);
26616 }
26617 computePyra(&splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0], M, &pyraVolume[6],
26618 pyraCentroid[6]);
26619 computePyra(&splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1], M, &pyraVolume[7],
26620 pyraCentroid[7]);
26621
26622
26623 computePoly5(p_1sss, p_2sss, p_2s, p_3ss, p_3sss, &area_c, coordinates_c, normalVec_c);
26624 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[8], pyraCentroid[8]);
26626 for(
MInt dim = 0; dim < 3; dim++) {
26627 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26628 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26629 }
26630
26631 computePoly6(p_1s, p_1ss, p_3s, p_3ss, p_2s, p_2ss, &area_c, coordinates_c, normalVec_c);
26632 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[9], pyraCentroid[9]);
26635 for(
MInt dim = 0; dim < 3; dim++) {
26636 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
26637 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
26638 }
26639
26640 for(
MInt i = 0; i < 10; i++) {
26641 volume_C += pyraVolume[i];
26642 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26643 }
26644 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26645
26647 for(
MInt dim = 0; dim < 3; dim++) {
26648 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26649 }
26650
26651 faceCut[face1] = true;
26652 faceCut[face2] = true;
26653 faceCut[face3] = true;
26654 faceCut[face4] = true;
26655 faceCut[face5] = true;
26656 faceCut[face6] = true;
26657
26658 } else {
26659
26661
26662
26663 computeTri(p_3, p_3sss, p_3ss, &faceVolume[face1], faceCentroid[face1], normal[face1]);
26664 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26665 faceVolume_0[face1] = faceVolume[face1];
26666 computeTri(p_2, p_2s, p_2sss, &faceVolume[face1], faceCentroid[face1], normal[face1]);
26667 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26668
26669 computeTri(p_4, p_1ss, p_3s, &splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0]);
26670 computeTri(p_0, p_3sss, p_1sss, &splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1]);
26671 splitCorner1 = cornersMCtoSOLVER[disambPointer[currentSubCase][4]];
26672 splitCorner2 = cornersMCtoSOLVER[disambPointer[currentSubCase][0]];
26673
26674 computeTri(p_5, p_2ss, p_1s, &splitFaceVolume2[0], splitFaceCentroid2[0], splitFaceNormal2[0]);
26675 computeTri(p_0, p_1sss, p_2sss, &splitFaceVolume2[1], splitFaceCentroid2[1], splitFaceNormal2[1]);
26676 splitCorner12 = cornersMCtoSOLVER[disambPointer[currentSubCase][5]];
26677 splitCorner22 = cornersMCtoSOLVER[disambPointer[currentSubCase][0]];
26678
26679 computeTri(p_1, p_1ss, p_1s, &faceVolume[face4], faceCentroid[face4], normal[face4]);
26680 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
26681 computeTri(p_2, p_2ss, p_2s, &faceVolume[face5], faceCentroid[face5], normal[face5]);
26682 correctFace(&faceVolume[face5], faceCentroid[face5], &faceVolume_0[face5], faceCentroid_0[face5]);
26683 computeTri(p_3, p_3ss, p_3s, &faceVolume[face6], faceCentroid[face6], normal[face6]);
26684 correctFace(&faceVolume[face6], faceCentroid[face6], &faceVolume_0[face6], faceCentroid_0[face6]);
26685
26686 splitFaceId = face2;
26687 splitFaceId2 = face3;
26688 maia::math::vecAvg<3>(M, p_1s, p_1ss, p_1sss, p_2s, p_2ss, p_2sss, p_3s, p_3ss, p_3sss);
26689
26690 faceVolume[face2] = F0;
26691 faceVolume[face3] = F0;
26692
26693 for(
MInt i = 0; i < 6; i++) {
26694 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]],
26695 M, &pyraVolume[i], pyraCentroid[i]);
26696 }
26697 computePyra(&splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0], M, &pyraVolume[6],
26698 pyraCentroid[6]);
26699 computePyra(&splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1], M, &pyraVolume[7],
26700 pyraCentroid[7]);
26701 computePyra(&splitFaceVolume2[0], splitFaceCentroid2[0], splitFaceNormal2[0], M, &pyraVolume[8],
26702 pyraCentroid[8]);
26703 computePyra(&splitFaceVolume2[1], splitFaceCentroid2[1], splitFaceNormal2[1], M, &pyraVolume[9],
26704 pyraCentroid[9]);
26705
26706 computePoly5(p_1ss, p_1sss, p_3sss, p_3ss, p_3s, &area_c, coordinates_c, normalVec_c);
26707 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[10], pyraCentroid[10]);
26709 for(
MInt dim = 0; dim < 3; dim++) {
26710 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26711 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26712 }
26713
26714 computePoly5(p_1sss, p_1s, p_2ss, p_2s, p_2sss, &area_c, coordinates_c, normalVec_c);
26715 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[11], pyraCentroid[11]);
26718 for(
MInt dim = 0; dim < 3; dim++) {
26719 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
26720 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
26721 }
26722
26723 computePoly3(p_1s, p_1sss, p_1ss, &area_c, coordinates_c, normalVec_c);
26724 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[12], pyraCentroid[12]);
26727 for(
MInt dim = 0; dim < 3; dim++) {
26728 m_bndryCells->
a[bndryId].m_srfcs[2]->m_coordinates[dim] = coordinates_c[dim];
26729 m_bndryCells->
a[bndryId].m_srfcs[2]->m_normalVector[dim] = normalVec_c[dim];
26730 }
26731
26732 for(
MInt i = 0; i < 13; i++) {
26733 volume_C += pyraVolume[i];
26734 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26735 }
26736 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26737
26739 for(
MInt dim = 0; dim < 3; dim++) {
26740 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26741 }
26742
26743 faceCut[face1] = true;
26744 faceCut[face2] = true;
26745 faceCut[face3] = true;
26746 faceCut[face4] = true;
26747 faceCut[face5] = true;
26748 faceCut[face6] = true;
26749 }
26750 break;
26751 }
26752 case 3:
26753 case 5:
26754 case 6: {
26755
26756
26757 if(disamb_tmp == 3) disambPointer_dummy = &tiling7_2_3STL[0][0];
26758 if(disamb_tmp == 5) disambPointer_dummy = &tiling7_2_5STL[0][0];
26759 if(disamb_tmp == 6) disambPointer_dummy = &tiling7_2_6STL[0][0];
26760
26761 for(
MInt i = 0; i < 16; i++)
26762 for(
MInt j = 0; j < 2; j++)
26763 disambPointer[i][j] = disambPointer_dummy[i * 2 + j];
26764
26765 for(
MInt i = 0; i < 6; i++) {
26766 faceVolume_0[i] = faceVolume[i];
26767 for(
MInt j = 0; j < 3; j++) {
26768 faceCentroid_0[i][j] = faceCentroid[i][j];
26769 }
26770 }
26771 cellVolume_0 = gridCellVolume;
26772 for(
MInt i = 0; i < 3; i++) {
26774 }
26775
26776 if(currentSubCase < 8) {
26777
26782 if(bndryId == 0)
26783 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! "
26784 << endl;
26787 } else {
26789 }
26790
26791 for(
MInt face = 0; face < 6; face++) {
26792 nfs_cur[face] = nfs7_1[currentSubCase][face];
26793 }
26794
26795
26796 subCaseDummy = disambPointer[currentSubCase][0];
26797 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
26798 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
26799 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26800 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
26801 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26802 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
26803 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26804 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
26805 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
26806 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
26807
26808 faceCut[face1] = true;
26809 faceCut[face2] = true;
26810 faceCut[face3] = true;
26811
26812 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
26813 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
26814 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
26815
26816 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
26817
26818 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
26819
26820 if(currentSubCase >= 8) {
26821 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26822 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26823 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26824
26825 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26826
26828 } else {
26829 for(
MInt face = 0; face < 6; face++) {
26830 nfs_cur[face] = nfs1[subCaseDummy][face];
26831 }
26832 }
26833
26834
26836 for(
MInt dim = 0; dim < 3; dim++) {
26837 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26838 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26839 }
26841 for(
MInt dim = 0; dim < 3; dim++) {
26842 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
26843 }
26844
26845 if(currentSubCase >= 8) {
26846 faceVolume_0[face1] = faceVolume[face1];
26847 faceVolume_0[face2] = faceVolume[face2];
26848 faceVolume_0[face3] = faceVolume[face3];
26849 cellVolume_0 = volume_C;
26850 volume_C = 0;
26851 for(
MInt i = 0; i < 3; i++) {
26852 cellCentroid_0[i] = coordinates_Cell[i];
26853 faceCentroid_0[face1][i] = faceCentroid[face1][i];
26854 faceCentroid_0[face2][i] = faceCentroid[face2][i];
26855 faceCentroid_0[face3][i] = faceCentroid[face3][i];
26856 coordinates_Cell[i] = 0;
26857 }
26858 } else {
26859 volume_C = 0;
26860 for(
MInt i = 0; i < 3; i++) {
26861 coordinates_Cell[i] = 0;
26862 }
26863 }
26864
26865
26866 subCaseDummy = disambPointer[currentSubCase][1];
26867
26868
26869 noFaces = 5;
26870
26871 p_0 = corner[cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][0]]];
26872 p_0s = corner[cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][1]]];
26873
26874 cutDummy = edgesMCtoSOLVER[tiling3_2STL[subCaseDummy][2]];
26875 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26876 cutDummy = edgesMCtoSOLVER[tiling3_2STL[subCaseDummy][3]];
26877 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26878 cutDummy = edgesMCtoSOLVER[tiling3_2STL[subCaseDummy][4]];
26879 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26880 cutDummy = edgesMCtoSOLVER[tiling3_2STL[subCaseDummy][5]];
26881 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26882 cutDummy = edgesMCtoSOLVER[tiling3_2STL[subCaseDummy][6]];
26883 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26884 cutDummy = edgesMCtoSOLVER[tiling3_2STL[subCaseDummy][7]];
26885 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
26886 face1 = facesMCtoSOLVER[tiling3_2STL[subCaseDummy][8]];
26887 face2 = facesMCtoSOLVER[tiling3_2STL[subCaseDummy][9]];
26888 face3 = facesMCtoSOLVER[tiling3_2STL[subCaseDummy][10]];
26889 face4 = facesMCtoSOLVER[tiling3_2STL[subCaseDummy][11]];
26890 face5 = facesMCtoSOLVER[tiling3_2STL[subCaseDummy][12]];
26891
26892 if(currentSubCase < 8) {
26893 for(
MInt face = 0; face < 6; face++) {
26894 nfs_cur_0[face] = nfs3_2[subCaseDummy][face];
26895 }
26896 faceCut_0[face1] = true;
26897 faceCut_0[face2] = true;
26898 faceCut_0[face3] = true;
26899 faceCut_0[face4] = true;
26900 faceCut_0[face5] = true;
26901
26902 computeTri(p_0, p_1, p_2, &faceVolume_0[face1], faceCentroid_0[face1], normal[face1]);
26903 computeTri(p_0, p_2, p_3, &faceVolume_0[face2], faceCentroid_0[face2], normal[face2]);
26904 computeTri(p_0s, p_2s, p_3s, &faceVolume_0[face3], faceCentroid_0[face3], normal[face3]);
26905 computeTri(p_0s, p_1s, p_2s, &faceVolume_0[face4], faceCentroid_0[face4], normal[face4]);
26906 computePoly6(p_0, p_1, p_3s, p_0s, p_1s, p_3, &faceVolume_0[face5], faceCentroid_0[face5],
26907 normal[face5]);
26908
26909 computePoly6(p_1, p_2, p_3, p_1s, p_2s, p_3s, &area_c, coordinates_c, normalVec_c);
26910
26911 maia::math::vecAvg<3>(M, p_0, p_0s, p_1, p_1s, p_2, p_2s, p_3, p_3s);
26912
26913 for(
MInt i = 0; i < 5; i++) {
26914 computePyra(&faceVolume_0[*facepointers[i]], faceCentroid_0[*facepointers[i]],
26915 normal[*facepointers[i]], M, &pyraVolume[i], pyraCentroid[i]);
26916 }
26917 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
26918
26919 for(
MInt i = 0; i < 6; i++) {
26920 volume_C += pyraVolume[i];
26921 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26922 }
26923 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26924
26925 absA = F0;
26926 for(
MInt i = 0; i < nDim; i++) {
26927 if(!nfs_cur_0[2 * i + 1]) faceVolume_0[2 * i + 1] = 0;
26928 if(!nfs_cur_0[2 * i]) faceVolume_0[2 * i] = 0;
26929 faceDiff[i] = faceVolume_0[2 * i + 1] - faceVolume_0[2 * i];
26930 absA +=
POW2(faceDiff[i]);
26931 }
26932
26933 for(
MInt i = 0; i < nDim; i++) {
26934 normalVec_c[i] = faceDiff[i] / sqrt(absA);
26935 }
26936 area_c = sqrt(absA);
26937
26938
26939
26941 for(
MInt dim = 0; dim < 3; dim++) {
26942 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
26943 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
26944 }
26945
26947 for(
MInt dim = 0; dim < 3; dim++) {
26948 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
26949 }
26950
26951
26952 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][0]]] = cellId2;
26953 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][1]]] = cellId2;
26954 } else {
26955 faceCut[face1] = true;
26956 faceCut[face2] = true;
26957 faceCut[face3] = true;
26958 faceCut[face4] = true;
26959 faceCut[face5] = true;
26960
26961 computeTri(p_0, p_1, p_2, &faceVolume[face1], faceCentroid[face1], normal[face1]);
26962 computeTri(p_0, p_2, p_3, &faceVolume[face2], faceCentroid[face2], normal[face2]);
26963 computeTri(p_0s, p_2s, p_3s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
26964 computeTri(p_0s, p_1s, p_2s, &faceVolume[face4], faceCentroid[face4], normal[face4]);
26965 computePoly6(p_0, p_1, p_3s, p_0s, p_1s, p_3, &faceVolume[face5], faceCentroid[face5], normal[face5]);
26966
26967 computePoly6(p_1, p_2, p_3, p_1s, p_2s, p_3s, &area_c, coordinates_c, normalVec_c);
26968
26969 maia::math::vecAvg<3>(M, p_0, p_0s, p_1, p_1s, p_2, p_2s, p_3, p_3s);
26970
26971 for(
MInt i = 0; i < 5; i++) {
26972 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]],
26973 M, &pyraVolume[i], pyraCentroid[i]);
26974 }
26975 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
26976
26977 for(
MInt i = 0; i < 6; i++) {
26978 volume_C += pyraVolume[i];
26979 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
26980 }
26981 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
26982
26983
26984 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
26985 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
26986 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
26987 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
26988 p_1ss = corner[cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][13]]];
26989 p_2ss = corner[cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][14]]];
26990 splitCorner1 = cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][13]];
26991 splitCorner2 = cornersMCtoSOLVER[tiling3_2STL[subCaseDummy][14]];
26992 computeTri(p_1s, p_1ss, p_3, &splitFaceVolume[0], splitFaceCentroid[0], splitFaceNormal[0]);
26993 computeTri(p_1, p_2ss, p_3s, &splitFaceVolume[1], splitFaceCentroid[1], splitFaceNormal[1]);
26994 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
26995
26997 splitFaceId = face5;
26998
26999
27002 for(
MInt dim = 0; dim < 3; dim++) {
27003 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
27004 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
27005 }
27006
27008 for(
MInt dim = 0; dim < 3; dim++) {
27009 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27010 }
27011 }
27012 break;
27013 }
27014 case 7: {
27015 if(currentSubCase < 8) {
27016
27025 if(bndryId == 0)
27026 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! "
27027 << endl;
27031 } else {
27033 }
27034
27035 for(
MInt face = 0; face < 6; face++) {
27036 nfs_cur[face] = nfs7_1[currentSubCase][face];
27037 }
27038
27039
27040 subCaseDummy = tiling7_1STL[currentSubCase][0];
27041 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
27042 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
27043 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27044 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
27045 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27046 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
27047 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27048 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
27049 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
27050 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
27051
27052 faceCut[face1] = true;
27053 faceCut[face2] = true;
27054 faceCut[face3] = true;
27055
27056 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
27057 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
27058 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
27059
27060 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27061
27062 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27063
27064 if(currentSubCase >= 8) {
27065 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27066 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27067 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27068
27069 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27070
27072 } else {
27073 for(
MInt face = 0; face < 6; face++) {
27074 nfs_cur[face] = nfs1[subCaseDummy][face];
27075 }
27076 }
27077
27078
27080 for(
MInt dim = 0; dim < 3; dim++) {
27081 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27082 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27083 }
27085 for(
MInt dim = 0; dim < 3; dim++) {
27086 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27087 }
27088
27089 if(currentSubCase >= 8) {
27090 faceVolume_0[face1] = faceVolume[face1];
27091 faceVolume_0[face2] = faceVolume[face2];
27092 faceVolume_0[face3] = faceVolume[face3];
27093 cellVolume_0 = volume_C;
27094 for(
MInt i = 0; i < 3; i++) {
27095 cellCentroid_0[i] = coordinates_Cell[i];
27096 faceCentroid_0[face1][i] = faceCentroid[face1][i];
27097 faceCentroid_0[face2][i] = faceCentroid[face2][i];
27098 faceCentroid_0[face3][i] = faceCentroid[face3][i];
27099 }
27100 }
27101
27102
27103 subCaseDummy = tiling7_1STL[currentSubCase][1];
27104 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
27105 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
27106 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27107 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
27108 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27109 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
27110 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27111 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
27112 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
27113 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
27114
27115 if(currentSubCase < 8) {
27116 for(
MInt face = 0; face < 6; face++) {
27117 nfs_cur_0[face] = nfs1[subCaseDummy][face];
27118 }
27119
27120 faceCut_0[face1] = true;
27121 faceCut_0[face2] = true;
27122 faceCut_0[face3] = true;
27123
27124 computeTri(p_0, p_1, p_3, &faceVolume_0[face1], faceCentroid_0[face1]);
27125 computeTri(p_0, p_3, p_2, &faceVolume_0[face2], faceCentroid_0[face2]);
27126 computeTri(p_0, p_2, p_1, &faceVolume_0[face3], faceCentroid_0[face3]);
27127
27128 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27129
27130 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27131
27132
27134 for(
MInt dim = 0; dim < 3; dim++) {
27135 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27136 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27137 }
27138
27140 for(
MInt dim = 0; dim < 3; dim++) {
27141 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
27142 }
27143
27144 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]] = cellId2;
27145 } else {
27146 faceCut[face1] = true;
27147 faceCut[face2] = true;
27148 faceCut[face3] = true;
27149
27150 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
27151 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
27152 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
27153
27154 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27155
27156 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27157
27158 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27159 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27160 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27161
27162 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27163
27165
27166
27169 for(
MInt dim = 0; dim < 3; dim++) {
27170 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
27171 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
27172 }
27173
27175 for(
MInt dim = 0; dim < 3; dim++) {
27176 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27177 }
27178 }
27179
27180 if(currentSubCase >= 8) {
27181 faceVolume_0[face1] = faceVolume[face1];
27182 faceVolume_0[face2] = faceVolume[face2];
27183 faceVolume_0[face3] = faceVolume[face3];
27184 cellVolume_0 = volume_C;
27185 for(
MInt i = 0; i < 3; i++) {
27186 cellCentroid_0[i] = coordinates_Cell[i];
27187 faceCentroid_0[face1][i] = faceCentroid[face1][i];
27188 faceCentroid_0[face2][i] = faceCentroid[face2][i];
27189 faceCentroid_0[face3][i] = faceCentroid[face3][i];
27190 }
27191 }
27192
27193
27194 subCaseDummy = tiling7_1STL[currentSubCase][2];
27195 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
27196 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
27197 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27198 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
27199 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27200 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
27201 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27202 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
27203 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
27204 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
27205
27206 if(currentSubCase < 8) {
27207 for(
MInt face = 0; face < 6; face++) {
27208 nfs_cur_00[face] = nfs1[subCaseDummy][face];
27209 }
27210
27211 faceCut_00[face1] = true;
27212 faceCut_00[face2] = true;
27213 faceCut_00[face3] = true;
27214
27215 computeTri(p_0, p_1, p_3, &faceVolume_00[face1], faceCentroid_00[face1]);
27216 computeTri(p_0, p_3, p_2, &faceVolume_00[face2], faceCentroid_00[face2]);
27217 computeTri(p_0, p_2, p_1, &faceVolume_00[face3], faceCentroid_00[face3]);
27218
27219 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27220
27221 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27222
27223
27225 for(
MInt dim = 0; dim < 3; dim++) {
27226 m_bndryCells->
a[bndryId3].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27227 m_bndryCells->
a[bndryId3].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27228 }
27229
27231 for(
MInt dim = 0; dim < 3; dim++) {
27232 m_bndryCells->
a[bndryId3].m_coordinates[dim] = coordinates_Cell[dim];
27233 }
27234
27235
27236 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]] = cellId3;
27237
27238 } else {
27239 faceCut[face1] = true;
27240 faceCut[face2] = true;
27241 faceCut[face3] = true;
27242
27243 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
27244 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
27245 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
27246
27247 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27248
27249 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27250
27251 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27252 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27253 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27254
27255 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27256
27258
27259
27262 for(
MInt dim = 0; dim < 3; dim++) {
27263 m_bndryCells->
a[bndryId].m_srfcs[2]->m_coordinates[dim] = coordinates_c[dim];
27264 m_bndryCells->
a[bndryId].m_srfcs[2]->m_normalVector[dim] = normalVec_c[dim];
27265 }
27266
27268 for(
MInt dim = 0; dim < 3; dim++) {
27269 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27270 }
27271 }
27272 break;
27273 }
27274 default: {
27275 stringstream errorMessage;
27276 errorMessage << "ERROR: Switch variable 'disamb_tmp' with value " << disamb_tmp
27277 << " not matching any case." << endl;
27278 mTerm(1, AT_, errorMessage.str());
27279 }
27280 }
27281 break;
27282 }
27283 case 10: {
27284 for(
MInt i = 0; i < 6; i++) {
27285 faceVolume_0[i] = faceVolume[i];
27286 for(
MInt j = 0; j < 3; j++) {
27287 faceCentroid_0[i][j] = faceCentroid[i][j];
27288 }
27289 }
27290 cellVolume_0 = gridCellVolume;
27291 for(
MInt i = 0; i < 3; i++) {
27293 }
27294
27295 if(disambiguation[bndryId] == 3) {
27296
27301 if(bndryId == 0)
27302 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! " << endl;
27305
27306
27307 subCaseDummy = tiling10_2STL[currentSubCase][0];
27308 for(
MInt face = 0; face < 6; face++) {
27309 nfs_cur[face] = nfs2[subCaseDummy][face];
27310 }
27311 p_0 = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][0]]];
27312 p_0s = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][1]]];
27313 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][2]];
27314 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27315 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][3]];
27316 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27317 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][4]];
27318 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27319 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][5]];
27320 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27321 face1 = facesMCtoSOLVER[tiling2STL[subCaseDummy][6]];
27322 face2 = facesMCtoSOLVER[tiling2STL[subCaseDummy][7]];
27323 face3 = facesMCtoSOLVER[tiling2STL[subCaseDummy][8]];
27324 face4 = facesMCtoSOLVER[tiling2STL[subCaseDummy][9]];
27325
27326 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2], normal[face2]);
27327 computeTri(p_0s, p_3s, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
27328 computeTrapez(p_0, p_0s, p_3s, p_3, &faceVolume[face1], faceCentroid[face1], normal[face1]);
27329 computeTrapez(p_2, p_2s, p_0s, p_0, &faceVolume[face4], faceCentroid[face4], normal[face4]);
27330
27331 computePoly4(p_3, p_3s, p_2s, p_2, &area_c, coordinates_c, normalVec_c);
27332
27333 maia::math::vecAvg<3>(M, p_0, p_0s, p_2, p_2s, p_3, p_3s);
27334
27335 for(
MInt i = 0; i < 4; i++) {
27336 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
27337 &pyraVolume[i], pyraCentroid[i]);
27338 }
27339 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[4], pyraCentroid[4]);
27340 for(
MInt i = 0; i < 5; i++) {
27341 volume_C += pyraVolume[i];
27342 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
27343 }
27344 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
27345
27346 faceCut[face1] = true;
27347 faceCut[face2] = true;
27348 faceCut[face3] = true;
27349 faceCut[face4] = true;
27350
27351
27353 for(
MInt dim = 0; dim < 3; dim++) {
27354 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27355 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27356 }
27357
27359 volume_C = 0;
27360 for(
MInt dim = 0; dim < 3; dim++) {
27361 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27362 coordinates_Cell[dim] = F0;
27363 }
27364
27365
27366 subCaseDummy = tiling10_2STL[currentSubCase][1];
27367 for(
MInt face = 0; face < 6; face++) {
27368 nfs_cur_0[face] = nfs2[subCaseDummy][face];
27369 }
27370 p_0 = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][0]]];
27371 p_0s = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][1]]];
27372 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][2]];
27373 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27374 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][3]];
27375 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27376 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][4]];
27377 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27378 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][5]];
27379 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27380 face1 = facesMCtoSOLVER[tiling2STL[subCaseDummy][6]];
27381 face2 = facesMCtoSOLVER[tiling2STL[subCaseDummy][7]];
27382 face3 = facesMCtoSOLVER[tiling2STL[subCaseDummy][8]];
27383 face4 = facesMCtoSOLVER[tiling2STL[subCaseDummy][9]];
27384
27385 computeTri(p_0, p_3, p_2, &faceVolume_0[face2], faceCentroid_0[face2], normal[face2]);
27386 computeTri(p_0s, p_3s, p_2s, &faceVolume_0[face3], faceCentroid_0[face3], normal[face3]);
27387 computeTrapez(p_0, p_0s, p_3s, p_3, &faceVolume_0[face1], faceCentroid_0[face1], normal[face1]);
27388 computeTrapez(p_2, p_2s, p_0s, p_0, &faceVolume_0[face4], faceCentroid_0[face4], normal[face4]);
27389
27390 computePoly4(p_3, p_3s, p_2s, p_2, &area_c, coordinates_c, normalVec_c);
27391
27392 maia::math::vecAvg<3>(M, p_0, p_0s, p_2, p_2s, p_3, p_3s);
27393
27394 for(
MInt i = 0; i < 4; i++) {
27395 computePyra(&faceVolume_0[*facepointers[i]], faceCentroid_0[*facepointers[i]], normal[*facepointers[i]],
27396 M, &pyraVolume[i], pyraCentroid[i]);
27397 }
27398 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[4], pyraCentroid[4]);
27399 for(
MInt i = 0; i < 5; i++) {
27400 volume_C += pyraVolume[i];
27401 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
27402 }
27403 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
27404
27405 faceCut_0[face1] = true;
27406 faceCut_0[face2] = true;
27407 faceCut_0[face3] = true;
27408 faceCut_0[face4] = true;
27409
27410
27412 for(
MInt dim = 0; dim < 3; dim++) {
27413 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27414 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27415 }
27416
27418 for(
MInt dim = 0; dim < 3; dim++) {
27419 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
27420 }
27421
27422
27423 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling2STL[subCaseDummy][0]]] = cellId2;
27424 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling2STL[subCaseDummy][1]]] = cellId2;
27425
27426 } else if(disambiguation[bndryId] == 0) {
27427
27428
27430 for(
MInt face = 0; face < 6; face++) {
27431 nfs_cur[face] = nfs10[currentSubCase][face];
27432 }
27433
27434
27435 subCaseDummy = tiling10_1STL[currentSubCase][0];
27436 p_0 = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][0]]];
27437 p_0s = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][1]]];
27438 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][2]];
27439 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27440 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][3]];
27441 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27442 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][4]];
27443 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27444 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][5]];
27445 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27446 face1 = facesMCtoSOLVER[tiling2STL[subCaseDummy][6]];
27447 face2 = facesMCtoSOLVER[tiling2STL[subCaseDummy][7]];
27448 face3 = facesMCtoSOLVER[tiling2STL[subCaseDummy][8]];
27449 face4 = facesMCtoSOLVER[tiling2STL[subCaseDummy][9]];
27450
27451 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2], normal[face2]);
27452 computeTri(p_0s, p_3s, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
27453 computeTrapez(p_0, p_0s, p_3s, p_3, &faceVolume[face1], faceCentroid[face1], normal[face1]);
27454 computeTrapez(p_2, p_2s, p_0s, p_0, &faceVolume[face4], faceCentroid[face4], normal[face4]);
27455
27456 computePoly4(p_3, p_3s, p_2s, p_2, &area_c, coordinates_c, normalVec_c);
27457
27458 maia::math::vecAvg<3>(M, p_0, p_0s, p_2, p_2s, p_3, p_3s);
27459
27460 for(
MInt i = 0; i < 4; i++) {
27461 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
27462 &pyraVolume[i], pyraCentroid[i]);
27463 }
27464 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[4], pyraCentroid[4]);
27465 for(
MInt i = 0; i < 5; i++) {
27466 volume_C += pyraVolume[i];
27467 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
27468 }
27469 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
27470
27471 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27472 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27473 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27474 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
27475
27476 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27477
27479
27480 faceCut[face1] = true;
27481 faceCut[face2] = true;
27482 faceCut[face3] = true;
27483 faceCut[face4] = true;
27484
27485
27487 for(
MInt dim = 0; dim < 3; dim++) {
27488 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27489 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27490 }
27491
27492 faceVolume_0[face1] = faceVolume[face1];
27493 faceVolume_0[face2] = faceVolume[face2];
27494 faceVolume_0[face3] = faceVolume[face3];
27495 faceVolume_0[face4] = faceVolume[face4];
27496 cellVolume_0 = volume_C;
27497 volume_C = 0;
27498 for(
MInt i = 0; i < 3; i++) {
27499 cellCentroid_0[i] = coordinates_Cell[i];
27500 faceCentroid_0[face1][i] = faceCentroid[face1][i];
27501 faceCentroid_0[face2][i] = faceCentroid[face2][i];
27502 faceCentroid_0[face3][i] = faceCentroid[face3][i];
27503 faceCentroid_0[face4][i] = faceCentroid[face4][i];
27504 coordinates_Cell[i] = 0;
27505 }
27506
27507
27508 subCaseDummy = tiling10_1STL[currentSubCase][1];
27509 p_0 = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][0]]];
27510 p_0s = corner[cornersMCtoSOLVER[tiling2STL[subCaseDummy][1]]];
27511 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][2]];
27512 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27513 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][3]];
27514 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27515 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][4]];
27516 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27517 cutDummy = edgesMCtoSOLVER[tiling2STL[subCaseDummy][5]];
27518 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27519 face1 = facesMCtoSOLVER[tiling2STL[subCaseDummy][6]];
27520 face2 = facesMCtoSOLVER[tiling2STL[subCaseDummy][7]];
27521 face3 = facesMCtoSOLVER[tiling2STL[subCaseDummy][8]];
27522 face4 = facesMCtoSOLVER[tiling2STL[subCaseDummy][9]];
27523
27524 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2], normal[face2]);
27525 computeTri(p_0s, p_3s, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
27526 computeTrapez(p_0, p_0s, p_3s, p_3, &faceVolume[face1], faceCentroid[face1], normal[face1]);
27527 computeTrapez(p_2, p_2s, p_0s, p_0, &faceVolume[face4], faceCentroid[face4], normal[face4]);
27528
27529 computePoly4(p_3, p_3s, p_2s, p_2, &area_c, coordinates_c, normalVec_c);
27530
27531 maia::math::vecAvg<3>(M, p_0, p_0s, p_2, p_2s, p_3, p_3s);
27532
27533 for(
MInt i = 0; i < 4; i++) {
27534 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
27535 &pyraVolume[i], pyraCentroid[i]);
27536 }
27537 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[4], pyraCentroid[4]);
27538 for(
MInt i = 0; i < 5; i++) {
27539 volume_C += pyraVolume[i];
27540 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
27541 }
27542 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
27543
27544 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27545 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27546 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27547 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
27548
27549 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27550
27552
27553 faceCut[face1] = true;
27554 faceCut[face2] = true;
27555 faceCut[face3] = true;
27556 faceCut[face4] = true;
27557
27558
27561 for(
MInt dim = 0; dim < 3; dim++) {
27562 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
27563 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
27564 }
27565
27567 for(
MInt dim = 0; dim < 3; dim++) {
27568 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27569 }
27570 }
27571 break;
27572 }
27573 case 12: {
27574 for(
MInt i = 0; i < 6; i++) {
27575 faceVolume_0[i] = faceVolume[i];
27576 for(
MInt j = 0; j < 3; j++) {
27577 faceCentroid_0[i][j] = faceCentroid[i][j];
27578 }
27579 }
27580 cellVolume_0 = gridCellVolume;
27581 for(
MInt i = 0; i < 3; i++) {
27583 }
27584
27585 if(disambiguation[bndryId] == 3) {
27586
27591 if(bndryId == 0)
27592 cerr << " problem in FvBndryCnd3D::createCutFaceMGC, assuming cell 0 is no split cell failed! " << endl;
27595
27596
27597 subCaseDummy = tiling12_1STL[23 - currentSubCase][0];
27598 for(
MInt face = 0; face < 6; face++) {
27599 nfs_cur[face] = nfs5[subCaseDummy][face];
27600 }
27601 p_0 = corner[cornersMCtoSOLVER[tiling5STL[subCaseDummy][0]]];
27602 p_1 = corner[cornersMCtoSOLVER[tiling5STL[subCaseDummy][1]]];
27603 p_2 = corner[cornersMCtoSOLVER[tiling5STL[subCaseDummy][2]]];
27604 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][3]];
27605 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27606 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][4]];
27607 p_0ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27608 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][5]];
27609 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27610 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][6]];
27611 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27612 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][7]];
27613 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27614 face1 = facesMCtoSOLVER[tiling5STL[subCaseDummy][8]];
27615 face2 = facesMCtoSOLVER[tiling5STL[subCaseDummy][9]];
27616 face3 = facesMCtoSOLVER[tiling5STL[subCaseDummy][10]];
27617 face4 = facesMCtoSOLVER[tiling5STL[subCaseDummy][11]];
27618 face5 = facesMCtoSOLVER[tiling5STL[subCaseDummy][12]];
27619
27620 computePoly5(p_0, p_0ss, p_2ss, p_2, p_1, &faceVolume[face1], faceCentroid[face1], normal[face1]);
27621 computeTri(p_0ss, p_0, p_0s, &faceVolume[face2], faceCentroid[face2], normal[face2]);
27622 computeTrapez(p_0, p_1, p_1s, p_0s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
27623 computeTrapez(p_1, p_1s, p_2s, p_2, &faceVolume[face4], faceCentroid[face4], normal[face4]);
27624 computeTri(p_2s, p_2, p_2ss, &faceVolume[face5], faceCentroid[face5], normal[face5]);
27625
27626 computePoly5(p_0s, p_1s, p_2s, p_2ss, p_0ss, &area_c, coordinates_c, normalVec_c);
27627
27628 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_0s, p_1s, p_2s, p_0ss, p_2ss);
27629
27630 for(
MInt i = 0; i < 5; i++) {
27631 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
27632 &pyraVolume[i], pyraCentroid[i]);
27633 }
27634 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
27635 for(
MInt i = 0; i < 6; i++) {
27636 volume_C += pyraVolume[i];
27637 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
27638 }
27639 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
27640
27641 faceCut[face1] = true;
27642 faceCut[face2] = true;
27643 faceCut[face3] = true;
27644 faceCut[face4] = true;
27645 faceCut[face5] = true;
27646
27647
27649 for(
MInt dim = 0; dim < 3; dim++) {
27650 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27651 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27652 }
27653
27655 volume_C = 0;
27656 for(
MInt dim = 0; dim < 3; dim++) {
27657 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27659 }
27660
27661
27662 subCaseDummy = tiling12_1STL[23 - currentSubCase][1];
27663 for(
MInt face = 0; face < 6; face++) {
27664 nfs_cur_0[face] = nfs1[subCaseDummy][face];
27665 }
27666 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
27667 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
27668 p_1 =
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27669 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
27670 p_2 =
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27671 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
27672 p_3 =
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27673 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
27674 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
27675 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
27676
27677 faceCut_0[face1] = true;
27678 faceCut_0[face2] = true;
27679 faceCut_0[face3] = true;
27680
27681 computeTri(p_0, p_1, p_3, &faceVolume_0[face1], faceCentroid_0[face1]);
27682 computeTri(p_0, p_3, p_2, &faceVolume_0[face2], faceCentroid_0[face2]);
27683 computeTri(p_0, p_2, p_1, &faceVolume_0[face3], faceCentroid_0[face3]);
27684
27685 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27686
27687 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27688
27690 for(
MInt dim = 0; dim < 3; dim++) {
27691 m_bndryCells->
a[bndryId2].m_coordinates[dim] = coordinates_Cell[dim];
27692 }
27693
27694
27696 for(
MInt dim = 0; dim < 3; dim++) {
27697 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27698 m_bndryCells->
a[bndryId2].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27699 }
27700
27701
27702 cornerCellMapping[ambIds[bndryId] * 8 + cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]] = cellId2;
27703
27704 } else if(disambiguation[bndryId] == 0) {
27705
27706
27708 for(
MInt face = 0; face < 6; face++) {
27709 nfs_cur[face] = nfs12_1[currentSubCase][face];
27710 }
27711
27712
27713 subCaseDummy = tiling12_1STL[currentSubCase][0];
27714
27715 p_0 = corner[cornersMCtoSOLVER[tiling5STL[subCaseDummy][0]]];
27716 p_1 = corner[cornersMCtoSOLVER[tiling5STL[subCaseDummy][1]]];
27717 p_2 = corner[cornersMCtoSOLVER[tiling5STL[subCaseDummy][2]]];
27718 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][3]];
27719 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27720 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][4]];
27721 p_0ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27722 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][5]];
27723 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27724 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][6]];
27725 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27726 cutDummy = edgesMCtoSOLVER[tiling5STL[subCaseDummy][7]];
27727 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27728 face1 = facesMCtoSOLVER[tiling5STL[subCaseDummy][8]];
27729 face2 = facesMCtoSOLVER[tiling5STL[subCaseDummy][9]];
27730 face3 = facesMCtoSOLVER[tiling5STL[subCaseDummy][10]];
27731 face4 = facesMCtoSOLVER[tiling5STL[subCaseDummy][11]];
27732 face5 = facesMCtoSOLVER[tiling5STL[subCaseDummy][12]];
27733
27734 computePoly5(p_0, p_0ss, p_2ss, p_2, p_1, &faceVolume[face1], faceCentroid[face1], normal[face1]);
27735 computeTri(p_0ss, p_0, p_0s, &faceVolume[face2], faceCentroid[face2], normal[face2]);
27736 computeTrapez(p_0, p_1, p_1s, p_0s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
27737 computeTrapez(p_1, p_1s, p_2s, p_2, &faceVolume[face4], faceCentroid[face4], normal[face4]);
27738 computeTri(p_2s, p_2, p_2ss, &faceVolume[face5], faceCentroid[face5], normal[face5]);
27739
27740 computePoly5(p_0s, p_1s, p_2s, p_2ss, p_0ss, &area_c, coordinates_c, normalVec_c);
27741
27742 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_0s, p_1s, p_2s, p_0ss, p_2ss);
27743
27744 for(
MInt i = 0; i < 5; i++) {
27745 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
27746 &pyraVolume[i], pyraCentroid[i]);
27747 }
27748 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
27749 for(
MInt i = 0; i < 6; i++) {
27750 volume_C += pyraVolume[i];
27751 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
27752 }
27753 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
27754
27755 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27756 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27757 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27758 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
27759 correctFace(&faceVolume[face5], faceCentroid[face5], &faceVolume_0[face5], faceCentroid_0[face5]);
27760
27761 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27762
27764
27765 faceCut[face1] = true;
27766 faceCut[face2] = true;
27767 faceCut[face3] = true;
27768 faceCut[face4] = true;
27769 faceCut[face5] = true;
27770
27771
27773 for(
MInt dim = 0; dim < 3; dim++) {
27774 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
27775 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
27776 }
27777
27778 faceVolume_0[face1] = faceVolume[face1];
27779 faceVolume_0[face2] = faceVolume[face2];
27780 faceVolume_0[face3] = faceVolume[face3];
27781 faceVolume_0[face4] = faceVolume[face4];
27782 faceVolume_0[face5] = faceVolume[face5];
27783 cellVolume_0 = volume_C;
27784 volume_C = 0;
27785 for(
MInt i = 0; i < 3; i++) {
27786 cellCentroid_0[i] = coordinates_Cell[i];
27787 faceCentroid_0[face1][i] = faceCentroid[face1][i];
27788 faceCentroid_0[face2][i] = faceCentroid[face2][i];
27789 faceCentroid_0[face3][i] = faceCentroid[face3][i];
27790 faceCentroid_0[face4][i] = faceCentroid[face4][i];
27791 faceCentroid_0[face5][i] = faceCentroid[face5][i];
27792 coordinates_Cell[i] = 0;
27793 }
27794
27795
27796 subCaseDummy = tiling12_1STL[currentSubCase][1];
27797 p_0 = corner[cornersMCtoSOLVER[tiling1STL[subCaseDummy][0]]];
27798 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][1]];
27799 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27800 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][2]];
27801 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27802 cutDummy = edgesMCtoSOLVER[tiling1STL[subCaseDummy][3]];
27803 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
27804 face1 = facesMCtoSOLVER[tiling1STL[subCaseDummy][4]];
27805 face2 = facesMCtoSOLVER[tiling1STL[subCaseDummy][5]];
27806 face3 = facesMCtoSOLVER[tiling1STL[subCaseDummy][6]];
27807
27808 faceCut[face1] = true;
27809 faceCut[face2] = true;
27810 faceCut[face3] = true;
27811
27812 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
27813 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
27814 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
27815
27816 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
27817
27818 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
27819
27820 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
27821 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
27822 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
27823
27824 correctCell(&volume_C, coordinates_Cell, &cellVolume_0, cellCentroid_0);
27825
27827
27828
27831 for(
MInt dim = 0; dim < 3; dim++) {
27832 m_bndryCells->
a[bndryId].m_srfcs[1]->m_coordinates[dim] = coordinates_c[dim];
27833 m_bndryCells->
a[bndryId].m_srfcs[1]->m_normalVector[dim] = normalVec_c[dim];
27834 }
27835
27837 for(
MInt dim = 0; dim < 3; dim++) {
27838 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
27839 }
27840 }
27841 break;
27842 }
27843 default: {
27844 mTerm(1, AT_,
"Unknown case probably in disamb_tmp");
27845 }
27846 }
27847
27848 for(
MInt i = 0; i < nDim; i++) {
27850 }
27851 if(splitCell[bndryId]) {
27852 for(
MInt i = 0; i < nDim; i++) {
27854 }
27855 }
27856 if(currentCase == 7 && disamb_tmp == 7) {
27857 for(
MInt i = 0; i < nDim; i++) {
27859 }
27860 }
27861
27862
27863
27864 if(!splitFace[bndryId]) {
27865 for(
MInt face = 0; face < 6; face++) {
27867 if(!faceCut[face]) continue;
27868 sideId = face % 2;
27869 spaceId = face / 2;
27870 if(nfs_cur[face]) {
27873 else {
27877 else
27878 continue;
27879 } else
27880 continue;
27881 }
27883 otherSideId = (sideId + 1) % 2;
27889 for(
MInt dim = 0; dim < nDim; dim++) {
27891 }
27894 } else {
27895 stringstream errorMessage;
27896 errorMessage << "Error in FvBndryCnd3D::createCutFaceMGC - 1. Problem with face " << face << endl;
27897 errorMessage << "bndryId: " << bndryId << endl;
27898 errorMessage <<
"cellId: " <<
cellId << endl;
27899 errorMessage <<
"bndryCnd: " <<
m_bndryCells->
a[bndryId].m_srfcs[0]->m_bndryCndId;
27900 errorMessage << "case, subcase: " << currentCase << ", " << currentSubCase << endl;
27901 errorMessage << "disamb: " << disambiguation[bndryId] << endl;
27902 errorMessage << "face Volumes: " << faceVolume[0] << ", " << faceVolume[1] << ", " << faceVolume[2] << ", "
27903 << faceVolume[3] << ", " << faceVolume[4] << ", " << faceVolume[5] << endl;
27904 errorMessage << "faceCentroids: ";
27905 for(
MInt i = 0; i < 6; i++) {
27906 errorMessage << faceCentroid[i][0] << ", " << faceCentroid[i][1] << ", " << faceCentroid[i][2] << "; "
27907 << endl;
27908 }
27909 errorMessage << "area_c = " << area_c << endl;
27910 errorMessage << "normalVec_c = " << normalVec_c[0] << ", " << normalVec_c[1] << ", " << normalVec_c[2]
27911 << endl;
27912 errorMessage << "coordinates_c = " << coordinates_c[0] << ", " << coordinates_c[1] << ", "
27913 << coordinates_c[2] << endl;
27914 errorMessage << "coordinates_Cell = " << coordinates_Cell[0] << ", " << coordinates_Cell[1] << ", "
27915 << coordinates_Cell[2] << endl;
27916 errorMessage << "volume_C = " << volume_C << endl;
27917 mTerm(1, AT_, errorMessage.str());
27918 }
27919 }
27920 }
27921 if(splitFace[bndryId]) {
27922 for(
MInt face = 0; face < 6; face++) {
27924 if(!faceCut[face]) continue;
27925 sideId = face % 2;
27926 spaceId = face / 2;
27927 if(nfs_cur[face]) {
27930 else {
27934 else
27935 continue;
27936 } else
27937 continue;
27938 }
27940 otherSideId = (sideId + 1) % 2;
27941 if(face == splitFaceId) {
27942
27948 for(
MInt dim = 0; dim < nDim; dim++) {
27950 }
27954 m_splitSurfaces[splitSurfaceIndexCounter * 6 + 2] = splitCorner1 + neighborCorner[face];
27955
27961 for(
MInt dim = 0; dim < nDim; dim++) {
27963 }
27966 m_splitSurfaces[splitSurfaceIndexCounter * 6 + 3 + 1] = splitCorner2;
27967 m_splitSurfaces[splitSurfaceIndexCounter * 6 + 3 + 2] = splitCorner2 + neighborCorner[face];
27968
27969 m_bndryCells->
a[bndryId].m_associatedSrfc[face] = -splitSurfaceIndexCounter;
27970 splitSurfaceIndexCounter++;
27971 } else {
27977 for(
MInt dim = 0; dim < nDim; dim++) {
27979 }
27982 }
27983 } else {
27984 stringstream errorMessage;
27985 errorMessage << "Error in FvBndryCnd3D::createCutFaceMGC - 2. Problem with face " << face << endl;
27986 errorMessage << "bndryId: " << bndryId << endl;
27987 errorMessage <<
"cellId: " <<
cellId << endl;
27988 errorMessage <<
"bndryCnd: " <<
m_bndryCells->
a[bndryId].m_srfcs[0]->m_bndryCndId;
27989 errorMessage << "case, subcase: " << currentCase << ", " << currentSubCase << endl;
27990 errorMessage << "disamb: " << disambiguation[bndryId] << endl;
27991 errorMessage << "face Volumes: " << faceVolume[0] << ", " << faceVolume[1] << ", " << faceVolume[2] << ", "
27992 << faceVolume[3] << ", " << faceVolume[4] << ", " << faceVolume[5] << endl;
27993 errorMessage << "faceCentroids: ";
27994 for(
MInt i = 0; i < 6; i++) {
27995 errorMessage << faceCentroid[i][0] << ", " << faceCentroid[i][1] << ", " << faceCentroid[i][2] << "; "
27996 << endl;
27997 }
27998 errorMessage << "area_c = " << area_c << endl;
27999 errorMessage << "normalVec_c = " << normalVec_c[0] << ", " << normalVec_c[1] << ", " << normalVec_c[2]
28000 << endl;
28001 errorMessage << "coordinates_c = " << coordinates_c[0] << ", " << coordinates_c[1] << ", "
28002 << coordinates_c[2] << endl;
28003 errorMessage << "coordinates_Cell = " << coordinates_Cell[0] << ", " << coordinates_Cell[1] << ", "
28004 << coordinates_Cell[2] << endl;
28005 errorMessage << "volume_C = " << volume_C << endl;
28006 mTerm(1, AT_, errorMessage.str());
28007 }
28008 }
28009 }
28010 if(currentCase == 7 && currentSubCase >= 8
28011 && (disamb_tmp == 1 || disamb_tmp == 2
28012 || disamb_tmp == 4)) {
28013
28014 for(
MInt face = 0; face < 6; face++) {
28015 if(face != splitFaceId2) continue;
28016 if(!faceCut[face]) continue;
28017 sideId = face % 2;
28018 spaceId = face / 2;
28019 if(nfs_cur[face]) {
28022 else {
28026 else
28027 continue;
28028 } else
28029 continue;
28030 }
28032 otherSideId = (sideId + 1) % 2;
28033
28039 for(
MInt dim = 0; dim < nDim; dim++) {
28041 }
28045 m_splitSurfaces[splitSurfaceIndexCounter * 6 + 2] = splitCorner12 + neighborCorner[face];
28046
28052 for(
MInt dim = 0; dim < nDim; dim++) {
28054 }
28057 m_splitSurfaces[splitSurfaceIndexCounter * 6 + 3 + 1] = splitCorner22;
28058 m_splitSurfaces[splitSurfaceIndexCounter * 6 + 3 + 2] = splitCorner22 + neighborCorner[face];
28059
28060 m_bndryCells->
a[bndryId].m_associatedSrfc[face] = -splitSurfaceIndexCounter;
28061 splitSurfaceIndexCounter++;
28062 } else {
28063 stringstream errorMessage;
28064 errorMessage << "Error in FvBndryCnd3D::createCutFaceMGC - 3. Problem with face " << face << endl;
28065 errorMessage << "bndryId: " << bndryId << endl;
28066 errorMessage <<
"cellId: " <<
cellId << endl;
28067 errorMessage <<
"bndryCnd: " <<
m_bndryCells->
a[bndryId].m_srfcs[0]->m_bndryCndId;
28068 errorMessage << "case, subcase: " << currentCase << ", " << currentSubCase << endl;
28069 errorMessage << "disamb: " << disambiguation[bndryId] << endl;
28070 errorMessage << "face Volumes: " << faceVolume[0] << ", " << faceVolume[1] << ", " << faceVolume[2] << ", "
28071 << faceVolume[3] << ", " << faceVolume[4] << ", " << faceVolume[5] << endl;
28072 errorMessage << "faceCentroids: ";
28073 for(
MInt i = 0; i < 6; i++) {
28074 errorMessage << faceCentroid[i][0] << ", " << faceCentroid[i][1] << ", " << faceCentroid[i][2] << "; "
28075 << endl;
28076 }
28077 errorMessage << "area_c = " << area_c << endl;
28078 errorMessage << "normalVec_c = " << normalVec_c[0] << ", " << normalVec_c[1] << ", " << normalVec_c[2]
28079 << endl;
28080 errorMessage << "coordinates_c = " << coordinates_c[0] << ", " << coordinates_c[1] << ", "
28081 << coordinates_c[2] << endl;
28082 errorMessage << "coordinates_Cell = " << coordinates_Cell[0] << ", " << coordinates_Cell[1] << ", "
28083 << coordinates_Cell[2] << endl;
28084 errorMessage << "volume_C = " << volume_C << endl;
28085 mTerm(1, AT_, errorMessage.str());
28086 }
28087 }
28088 }
28089 if(splitCell[bndryId]) {
28090 for(
MInt face = 0; face < 6; face++) {
28092 if(!faceCut_0[face]) continue;
28093 sideId = face % 2;
28094 spaceId = face / 2;
28095 if(nfs_cur_0[face]) {
28096
28097 MInt splitChildId = cellId2;
28102 else {
28106 else
28107 continue;
28108 } else
28109 continue;
28110 }
28111 cellId2 = splitChildId;
28113 otherSideId = (sideId + 1) % 2;
28119 for(
MInt dim = 0; dim < nDim; dim++) {
28121 }
28124 } else {
28125 stringstream errorMessage;
28126 errorMessage << "Error in FvBndryCnd3D::createCutFaceMGC - 4. Problem with face " << face << endl;
28127 errorMessage << "bndryId: " << bndryId2 << endl;
28128 errorMessage << "cellId: " << cellId2 << endl;
28129 errorMessage << "split cell! Twin: " << bndryId << endl;
28130 errorMessage <<
"bndryCnd: " <<
m_bndryCells->
a[bndryId2].m_srfcs[0]->m_bndryCndId;
28131 errorMessage << "case, subcase: " << currentCase << ", " << currentSubCase << endl;
28132 errorMessage << "disamb: " << disambiguation[bndryId] << endl;
28133 errorMessage << "face Volumes: " << faceVolume_0[0] << ", " << faceVolume_0[1] << ", " << faceVolume_0[2]
28134 << ", " << faceVolume_0[3] << ", " << faceVolume_0[4] << ", " << faceVolume_0[5] << endl;
28135 errorMessage << "faceCentroids: ";
28136 for(
MInt i = 0; i < 6; i++) {
28137 errorMessage << faceCentroid_0[i][0] << ", " << faceCentroid_0[i][1] << ", " << faceCentroid_0[i][2]
28138 << "; " << endl;
28139 }
28140 errorMessage << "area_c = " << area_c << endl;
28141 errorMessage << "normalVec_c = " << normalVec_c[0] << ", " << normalVec_c[1] << ", " << normalVec_c[2]
28142 << endl;
28143 errorMessage << "coordinates_c = " << coordinates_c[0] << ", " << coordinates_c[1] << ", "
28144 << coordinates_c[2] << endl;
28145 errorMessage << "coordinates_Cell = " << coordinates_Cell[0] << ", " << coordinates_Cell[1] << ", "
28146 << coordinates_Cell[2] << endl;
28147 errorMessage << "volume_C = " << volume_C << endl;
28148 mTerm(1, AT_, errorMessage.str());
28149 }
28150 }
28151 }
28152
28153 if(currentCase == 7 && disamb_tmp == 7) {
28154 for(
MInt face = 0; face < 6; face++) {
28156 if(!faceCut_00[face]) continue;
28157 sideId = face % 2;
28158 spaceId = face / 2;
28159 if(nfs_cur_00[face]) {
28160
28161 MInt splitChildId = cellId3;
28166 else {
28170 else
28171 continue;
28172 } else
28173 continue;
28174 }
28175 cellId3 = splitChildId;
28177 otherSideId = (sideId + 1) % 2;
28183 for(
MInt dim = 0; dim < nDim; dim++) {
28185 }
28188 } else {
28189 stringstream errorMessage;
28190 errorMessage << "Error in FvBndryCnd3D::createCutFaceMGC - 5. Problem with face " << face << endl;
28191 errorMessage << "bndryId: " << bndryId3 << endl;
28192 errorMessage << "cellId: " << cellId3 << endl;
28193 errorMessage << "split cell! Twin (triple): " << bndryId << endl;
28194 errorMessage <<
"bndryCnd: " <<
m_bndryCells->
a[bndryId3].m_srfcs[0]->m_bndryCndId;
28195 errorMessage << "case, subcase: " << currentCase << ", " << currentSubCase << endl;
28196 errorMessage << "disamb: " << disambiguation[bndryId] << endl;
28197 errorMessage << "face Volumes: " << faceVolume_00[0] << ", " << faceVolume_00[1] << ", " << faceVolume_00[2]
28198 << ", " << faceVolume_00[3] << ", " << faceVolume_00[4] << ", " << faceVolume_00[5] << endl;
28199 errorMessage << "faceCentroids: ";
28200 for(
MInt i = 0; i < 6; i++) {
28201 errorMessage << faceCentroid_00[i][0] << ", " << faceCentroid_00[i][1] << ", " << faceCentroid_00[i][2]
28202 << "; " << endl;
28203 }
28204 errorMessage << "area_c = " << area_c << endl;
28205 errorMessage << "normalVec_c = " << normalVec_c[0] << ", " << normalVec_c[1] << ", " << normalVec_c[2]
28206 << endl;
28207 errorMessage << "coordinates_c = " << coordinates_c[0] << ", " << coordinates_c[1] << ", "
28208 << coordinates_c[2] << endl;
28209 errorMessage << "coordinates_Cell = " << coordinates_Cell[0] << ", " << coordinates_Cell[1] << ", "
28210 << coordinates_Cell[2] << endl;
28211 mTerm(1, AT_, errorMessage.str());
28212 }
28213 }
28214 }
28215
28216 for(
MInt face = 0; face < 6; face++) {
28217 if(!nfs_cur[face]) {
28219 }
28220 }
28221 if(splitCell[bndryId]) {
28222 for(
MInt face = 0; face < 6; face++) {
28223 if(!nfs_cur_0[face]) {
28225 }
28226 }
28227 }
28228 if(currentCase == 7 && disamb_tmp == 7) {
28229 for(
MInt face = 0; face < 6; face++) {
28230 if(!nfs_cur_00[face]) {
28232 }
28233 }
28234 }
28235 } else {
28237 stringstream fileName;
28238 switch(currentCase) {
28239 case 0:
28240 cerr << "FvBndryCnd3D::createCutFaceMGC - Error: Cell is not a boundary cell!" << endl;
28241 cerr <<
"CellId: " <<
cellId <<
" BndryId: " << bndryId << endl;
28243 fileName <<
cellId <<
".stl";
28245 cerr << "Wrote stl-file of cell. FileName: " << (fileName.str()) << endl;
28246 break;
28247 case 1: {
28248
28249 for(
MInt face = 0; face < 6; face++) {
28250 nfs_cur[face] = nfs1[currentSubCase][face];
28251 }
28252 p_0 = corner[cornersMCtoSOLVER[tiling1STL[currentSubCase][0]]];
28253 noFaces = 3;
28254 if(currentSubCase < 8) {
28255 } else {
28256 for(
MInt i = 0; i < 6; i++) {
28257 faceVolume_0[i] = faceVolume[i];
28258 for(
MInt j = 0; j < 3; j++) {
28259 faceCentroid_0[i][j] = faceCentroid[i][j];
28260 }
28261 }
28262 }
28263 cutDummy = edgesMCtoSOLVER[tiling1STL[currentSubCase][1]];
28264 p_1 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28265 cutDummy = edgesMCtoSOLVER[tiling1STL[currentSubCase][2]];
28266 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28267 cutDummy = edgesMCtoSOLVER[tiling1STL[currentSubCase][3]];
28268 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28269 face1 = facesMCtoSOLVER[tiling1STL[currentSubCase][4]];
28270 face2 = facesMCtoSOLVER[tiling1STL[currentSubCase][5]];
28271 face3 = facesMCtoSOLVER[tiling1STL[currentSubCase][6]];
28272
28273 computeTri(p_0, p_1, p_3, &faceVolume[face1], faceCentroid[face1]);
28274 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2]);
28275 computeTri(p_0, p_2, p_1, &faceVolume[face3], faceCentroid[face3]);
28276
28277 computePoly3(p_1, p_2, p_3, &area_c, coordinates_c, normalVec_c);
28278
28279 computeTetra(p_0, p_1, p_2, p_3, &volume_C, coordinates_Cell);
28280
28281 if(currentSubCase >= 8) {
28282 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
28283 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
28284 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
28285
28287
28289 }
28290
28292
28294 for(
MInt dim = 0; dim < 3; dim++) {
28295 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28296 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28297 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28298 }
28299 } break;
28300 case 2: {
28301
28302 for(
MInt face = 0; face < 6; face++) {
28303 nfs_cur[face] = nfs2[currentSubCase][face];
28304 }
28305 noFaces = 4;
28306 p_0 = corner[cornersMCtoSOLVER[tiling2STL[currentSubCase][0]]];
28307 p_0s = corner[cornersMCtoSOLVER[tiling2STL[currentSubCase][1]]];
28308 if(currentSubCase < 12) {
28309 } else {
28310 for(
MInt i = 0; i < 6; i++) {
28311 faceVolume_0[i] = faceVolume[i];
28312 for(
MInt j = 0; j < 3; j++) {
28313 faceCentroid_0[i][j] = faceCentroid[i][j];
28314 }
28315 }
28316 }
28317
28318 cutDummy = edgesMCtoSOLVER[tiling2STL[currentSubCase][2]];
28319 p_3 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28320 cutDummy = edgesMCtoSOLVER[tiling2STL[currentSubCase][3]];
28321 p_2 =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28322 cutDummy = edgesMCtoSOLVER[tiling2STL[currentSubCase][4]];
28323 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28324 cutDummy = edgesMCtoSOLVER[tiling2STL[currentSubCase][5]];
28325 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28326 face1 = facesMCtoSOLVER[tiling2STL[currentSubCase][6]];
28327 face2 = facesMCtoSOLVER[tiling2STL[currentSubCase][7]];
28328 face3 = facesMCtoSOLVER[tiling2STL[currentSubCase][8]];
28329 face4 = facesMCtoSOLVER[tiling2STL[currentSubCase][9]];
28330
28331 computeTri(p_0, p_3, p_2, &faceVolume[face2], faceCentroid[face2], normal[face2]);
28332 computeTri(p_0s, p_3s, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
28333 computeTrapez(p_0, p_0s, p_3s, p_3, &faceVolume[face1], faceCentroid[face1], normal[face1]);
28334 computeTrapez(p_2, p_2s, p_0s, p_0, &faceVolume[face4], faceCentroid[face4], normal[face4]);
28335
28336 computePoly4(p_3, p_3s, p_2s, p_2, &area_c, coordinates_c, normalVec_c);
28337
28338 maia::math::vecAvg<3>(M, p_0, p_0s, p_2, p_2s, p_3, p_3s);
28339
28340 for(
MInt i = 0; i < 4; i++) {
28341 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
28342 &pyraVolume[i], pyraCentroid[i]);
28343 }
28344 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[4], pyraCentroid[4]);
28345 for(
MInt i = 0; i < 5; i++) {
28346 volume_C += pyraVolume[i];
28347 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
28348 }
28349 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
28350
28351 if(currentSubCase >= 12) {
28352 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
28353 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
28354 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
28355 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
28356
28358
28360 }
28361
28363
28365 for(
MInt dim = 0; dim < 3; dim++) {
28366 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28367 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28368 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28369 }
28370 } break;
28371 case 5: {
28372
28373 for(
MInt face = 0; face < 6; face++) {
28374 nfs_cur[face] = nfs5[currentSubCase][face];
28375 }
28376 noFaces = 5;
28377 p_0 = corner[cornersMCtoSOLVER[tiling5STL[currentSubCase][0]]];
28378 p_1 = corner[cornersMCtoSOLVER[tiling5STL[currentSubCase][1]]];
28379 p_2 = corner[cornersMCtoSOLVER[tiling5STL[currentSubCase][2]]];
28380 if(currentSubCase < 24) {
28381 } else {
28382 for(
MInt i = 0; i < 6; i++) {
28383 faceVolume_0[i] = faceVolume[i];
28384 for(
MInt j = 0; j < 3; j++) {
28385 faceCentroid_0[i][j] = faceCentroid[i][j];
28386 }
28387 }
28388 }
28389
28390 cutDummy = edgesMCtoSOLVER[tiling5STL[currentSubCase][3]];
28391 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28392 cutDummy = edgesMCtoSOLVER[tiling5STL[currentSubCase][4]];
28393 p_0ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28394 cutDummy = edgesMCtoSOLVER[tiling5STL[currentSubCase][5]];
28395 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28396 cutDummy = edgesMCtoSOLVER[tiling5STL[currentSubCase][6]];
28397 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28398 cutDummy = edgesMCtoSOLVER[tiling5STL[currentSubCase][7]];
28399 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28400 face1 = facesMCtoSOLVER[tiling5STL[currentSubCase][8]];
28401 face2 = facesMCtoSOLVER[tiling5STL[currentSubCase][9]];
28402 face3 = facesMCtoSOLVER[tiling5STL[currentSubCase][10]];
28403 face4 = facesMCtoSOLVER[tiling5STL[currentSubCase][11]];
28404 face5 = facesMCtoSOLVER[tiling5STL[currentSubCase][12]];
28405
28406 computePoly5(p_0, p_0ss, p_2ss, p_2, p_1, &faceVolume[face1], faceCentroid[face1], normal[face1]);
28407 computeTri(p_0ss, p_0, p_0s, &faceVolume[face2], faceCentroid[face2], normal[face2]);
28408 computeTrapez(p_0, p_1, p_1s, p_0s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
28409 computeTrapez(p_1, p_1s, p_2s, p_2, &faceVolume[face4], faceCentroid[face4], normal[face4]);
28410 computeTri(p_2s, p_2, p_2ss, &faceVolume[face5], faceCentroid[face5], normal[face5]);
28411
28412 computePoly5(p_0s, p_1s, p_2s, p_2ss, p_0ss, &area_c, coordinates_c, normalVec_c);
28413
28414 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_0s, p_1s, p_2s, p_0ss, p_2ss);
28415
28416 for(
MInt i = 0; i < 5; i++) {
28417 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
28418 &pyraVolume[i], pyraCentroid[i]);
28419 }
28420 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
28421 for(
MInt i = 0; i < 6; i++) {
28422 volume_C += pyraVolume[i];
28423 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
28424 }
28425 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
28426
28427 if(currentSubCase >= 24) {
28428 correctFace(&faceVolume[face1], faceCentroid[face1], &faceVolume_0[face1], faceCentroid_0[face1]);
28429 correctFace(&faceVolume[face2], faceCentroid[face2], &faceVolume_0[face2], faceCentroid_0[face2]);
28430 correctFace(&faceVolume[face3], faceCentroid[face3], &faceVolume_0[face3], faceCentroid_0[face3]);
28431 correctFace(&faceVolume[face4], faceCentroid[face4], &faceVolume_0[face4], faceCentroid_0[face4]);
28432 correctFace(&faceVolume[face5], faceCentroid[face5], &faceVolume_0[face5], faceCentroid_0[face5]);
28433
28435
28437 }
28438
28440
28442 for(
MInt dim = 0; dim < 3; dim++) {
28443 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28444 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28445 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28446 }
28447 } break;
28448 case 8: {
28449
28450 for(
MInt face = 0; face < 6; face++) {
28451 nfs_cur[face] = nfs8[currentSubCase][face];
28452 }
28453 noFaces = 4;
28454 p_0 = corner[cornersMCtoSOLVER[tiling8STL[currentSubCase][0]]];
28455 p_1 = corner[cornersMCtoSOLVER[tiling8STL[currentSubCase][1]]];
28456 p_2 = corner[cornersMCtoSOLVER[tiling8STL[currentSubCase][2]]];
28457 p_3 = corner[cornersMCtoSOLVER[tiling8STL[currentSubCase][3]]];
28458
28459 cutDummy = edgesMCtoSOLVER[tiling8STL[currentSubCase][4]];
28460 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28461 cutDummy = edgesMCtoSOLVER[tiling8STL[currentSubCase][5]];
28462 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28463 cutDummy = edgesMCtoSOLVER[tiling8STL[currentSubCase][6]];
28464 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28465 cutDummy = edgesMCtoSOLVER[tiling8STL[currentSubCase][7]];
28466 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28467 face5 = facesMCtoSOLVER[tiling8STL[currentSubCase][8]];
28468 face1 = facesMCtoSOLVER[tiling8STL[currentSubCase][9]];
28469 face2 = facesMCtoSOLVER[tiling8STL[currentSubCase][10]];
28470 face3 = facesMCtoSOLVER[tiling8STL[currentSubCase][11]];
28471 face4 = facesMCtoSOLVER[tiling8STL[currentSubCase][12]];
28472
28473 computeTrapez(p_0, p_1, p_1s, p_0s, &faceVolume[face1], faceCentroid[face1], normal[face1]);
28474 computeTrapez(p_1, p_2, p_2s, p_1s, &faceVolume[face2], faceCentroid[face2], normal[face2]);
28475 computeTrapez(p_2, p_3, p_3s, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
28476 computeTrapez(p_3, p_0, p_0s, p_3s, &faceVolume[face4], faceCentroid[face4], normal[face4]);
28477
28478 computeTrapez(p_0, p_3, p_2, p_1, &faceVolume[face5], faceCentroid[face5], normal[face5]);
28479
28480 computePoly4(p_0s, p_1s, p_2s, p_3s, &area_c, coordinates_c, normalVec_c);
28481
28482 maia::math::vecAvg<3>(M, p_0, p_0s, p_1, p_1s, p_2, p_2s, p_3, p_3s);
28483
28484 for(
MInt i = 0; i < 5; i++) {
28485 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
28486 &pyraVolume[i], pyraCentroid[i]);
28487 }
28488 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[5], pyraCentroid[5]);
28489 for(
MInt i = 0; i < 6; i++) {
28490 volume_C += pyraVolume[i];
28491 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
28492 }
28493 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
28494
28496
28498 for(
MInt dim = 0; dim < 3; dim++) {
28499 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28500 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28501 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28502 }
28503 } break;
28504 case 9: {
28505
28506 for(
MInt face = 0; face < 6; face++) {
28507 nfs_cur[face] = nfs9[currentSubCase][face];
28508 }
28509 noFaces = 6;
28510 p_0 = corner[cornersMCtoSOLVER[tiling9STL[currentSubCase][0]]];
28511 p_1 = corner[cornersMCtoSOLVER[tiling9STL[currentSubCase][1]]];
28512 p_2 = corner[cornersMCtoSOLVER[tiling9STL[currentSubCase][2]]];
28513 p_3 = corner[cornersMCtoSOLVER[tiling9STL[currentSubCase][3]]];
28514
28515 cutDummy = edgesMCtoSOLVER[tiling9STL[currentSubCase][4]];
28516 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28517 cutDummy = edgesMCtoSOLVER[tiling9STL[currentSubCase][5]];
28518 p_1ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28519 cutDummy = edgesMCtoSOLVER[tiling9STL[currentSubCase][6]];
28520 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28521 cutDummy = edgesMCtoSOLVER[tiling9STL[currentSubCase][7]];
28522 p_2ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28523 cutDummy = edgesMCtoSOLVER[tiling9STL[currentSubCase][8]];
28524 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28525 cutDummy = edgesMCtoSOLVER[tiling9STL[currentSubCase][9]];
28526 p_3ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28527 face1 = facesMCtoSOLVER[tiling9STL[currentSubCase][10]];
28528 face2 = facesMCtoSOLVER[tiling9STL[currentSubCase][11]];
28529 face3 = facesMCtoSOLVER[tiling9STL[currentSubCase][12]];
28530 face4 = facesMCtoSOLVER[tiling9STL[currentSubCase][13]];
28531 face5 = facesMCtoSOLVER[tiling9STL[currentSubCase][14]];
28532 face6 = facesMCtoSOLVER[tiling9STL[currentSubCase][15]];
28533
28534 computePoly5(p_0, p_3, p_3ss, p_2s, p_2, &faceVolume[face1], faceCentroid[face1], normal[face1]);
28535 computePoly5(p_0, p_1, p_1ss, p_3s, p_3, &faceVolume[face3], faceCentroid[face3], normal[face3]);
28536 computePoly5(p_0, p_2, p_2ss, p_1s, p_1, &faceVolume[face5], faceCentroid[face5], normal[face5]);
28537 computeTri(p_1, p_1s, p_1ss, &faceVolume[face2], faceCentroid[face2], normal[face2]);
28538 computeTri(p_2, p_2s, p_2ss, &faceVolume[face4], faceCentroid[face4], normal[face4]);
28539 computeTri(p_3, p_3s, p_3ss, &faceVolume[face6], faceCentroid[face6], normal[face6]);
28540
28541 computePoly6(p_1ss, p_1s, p_2ss, p_2s, p_3ss, p_3s, &area_c, coordinates_c, normalVec_c);
28542
28543 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_3, p_1s, p_2s, p_3s, p_1ss, p_2ss, p_3ss);
28544
28545 for(
MInt i = 0; i < 6; i++) {
28546 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
28547 &pyraVolume[i], pyraCentroid[i]);
28548 }
28549 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[6], pyraCentroid[6]);
28550 for(
MInt i = 0; i < 7; i++) {
28551 volume_C += pyraVolume[i];
28552 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
28553 }
28554 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
28555
28557
28559 for(
MInt dim = 0; dim < 3; dim++) {
28560 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28561 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28562 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28563 }
28564 }
28565
28566 break;
28567 case 11: {
28568
28569 for(
MInt face = 0; face < 6; face++) {
28570 nfs_cur[face] = nfs11[currentSubCase][face];
28571 }
28572 noFaces = 6;
28573 p_0 = corner[cornersMCtoSOLVER[tiling11STL[currentSubCase][0]]];
28574 p_1 = corner[cornersMCtoSOLVER[tiling11STL[currentSubCase][1]]];
28575 p_2 = corner[cornersMCtoSOLVER[tiling11STL[currentSubCase][2]]];
28576 p_3 = corner[cornersMCtoSOLVER[tiling11STL[currentSubCase][3]]];
28577
28578 cutDummy = edgesMCtoSOLVER[tiling11STL[currentSubCase][4]];
28579 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28580 cutDummy = edgesMCtoSOLVER[tiling11STL[currentSubCase][5]];
28581 p_0ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28582 cutDummy = edgesMCtoSOLVER[tiling11STL[currentSubCase][6]];
28583 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28584 cutDummy = edgesMCtoSOLVER[tiling11STL[currentSubCase][7]];
28585 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28586 cutDummy = edgesMCtoSOLVER[tiling11STL[currentSubCase][8]];
28587 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28588 cutDummy = edgesMCtoSOLVER[tiling11STL[currentSubCase][9]];
28589 p_3ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28590 face1 = facesMCtoSOLVER[tiling11STL[currentSubCase][10]];
28591 face2 = facesMCtoSOLVER[tiling11STL[currentSubCase][11]];
28592 face3 = facesMCtoSOLVER[tiling11STL[currentSubCase][12]];
28593 face4 = facesMCtoSOLVER[tiling11STL[currentSubCase][13]];
28594 face5 = facesMCtoSOLVER[tiling11STL[currentSubCase][14]];
28595 face6 = facesMCtoSOLVER[tiling11STL[currentSubCase][15]];
28596
28597 computeTri(p_0, p_0s, p_0ss, &faceVolume[face1], faceCentroid[face1], normal[face1]);
28598 computePoly5(p_3s, p_3, p_2, p_1, p_1s, &faceVolume[face2], faceCentroid[face2], normal[face2]);
28599 computeTrapez(p_2, p_3, p_3ss, p_2s, &faceVolume[face3], faceCentroid[face3], normal[face3]);
28600 computeTrapez(p_1, p_0, p_0ss, p_1s, &faceVolume[face4], faceCentroid[face4], normal[face4]);
28601 computePoly5(p_1, p_2, p_2s, p_0s, p_0, &faceVolume[face5], faceCentroid[face5], normal[face5]);
28602 computeTri(p_3, p_3s, p_3ss, &faceVolume[face6], faceCentroid[face6], normal[face6]);
28603
28604 computePoly6(p_0ss, p_0s, p_2s, p_3ss, p_3s, p_1s, &area_c, coordinates_c, normalVec_c);
28605
28606 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_3, p_0s, p_1s, p_2s, p_3s, p_0ss, p_3ss);
28607
28608 for(
MInt i = 0; i < 6; i++) {
28609 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
28610 &pyraVolume[i], pyraCentroid[i]);
28611 }
28612 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[6], pyraCentroid[6]);
28613 for(
MInt i = 0; i < 7; i++) {
28614 volume_C += pyraVolume[i];
28615 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
28616 }
28617 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
28618
28620
28622 for(
MInt dim = 0; dim < 3; dim++) {
28623 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28624 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28625 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28626 }
28627 break;
28628 }
28629 case 14: {
28630
28631 for(
MInt face = 0; face < 6; face++) {
28632 nfs_cur[face] = nfs14[currentSubCase][face];
28633 }
28634 noFaces = 6;
28635 p_0 = corner[cornersMCtoSOLVER[tiling14STL[currentSubCase][0]]];
28636 p_1 = corner[cornersMCtoSOLVER[tiling14STL[currentSubCase][1]]];
28637 p_2 = corner[cornersMCtoSOLVER[tiling14STL[currentSubCase][2]]];
28638 p_3 = corner[cornersMCtoSOLVER[tiling14STL[currentSubCase][3]]];
28639
28640 cutDummy = edgesMCtoSOLVER[tiling14STL[currentSubCase][4]];
28641 p_0s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28642 cutDummy = edgesMCtoSOLVER[tiling14STL[currentSubCase][5]];
28643 p_0ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28644 cutDummy = edgesMCtoSOLVER[tiling14STL[currentSubCase][6]];
28645 p_1s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28646 cutDummy = edgesMCtoSOLVER[tiling14STL[currentSubCase][7]];
28647 p_2s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28648 cutDummy = edgesMCtoSOLVER[tiling14STL[currentSubCase][8]];
28649 p_3s =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28650 cutDummy = edgesMCtoSOLVER[tiling14STL[currentSubCase][9]];
28651 p_3ss =
m_bndryCells->
a[bndryId].m_srfcs[0]->m_cutCoordinates[cutPoints[cutDummy]];
28652 face1 = facesMCtoSOLVER[tiling14STL[currentSubCase][10]];
28653 face2 = facesMCtoSOLVER[tiling14STL[currentSubCase][11]];
28654 face3 = facesMCtoSOLVER[tiling14STL[currentSubCase][12]];
28655 face4 = facesMCtoSOLVER[tiling14STL[currentSubCase][13]];
28656 face5 = facesMCtoSOLVER[tiling14STL[currentSubCase][14]];
28657 face6 = facesMCtoSOLVER[tiling14STL[currentSubCase][15]];
28658
28659 computeTri(p_0, p_0s, p_0ss, &faceVolume[face1], faceCentroid[face1], normal[face1]);
28660 computePoly5(p_2, p_3, p_3ss, p_1s, p_1, &faceVolume[face2], faceCentroid[face2], normal[face2]);
28661 computeTrapez(p_2s, p_3s, p_3, p_2, &faceVolume[face3], faceCentroid[face3], normal[face3]);
28662 computeTrapez(p_1s, p_0s, p_0, p_1, &faceVolume[face4], faceCentroid[face4], normal[face4]);
28663 computePoly5(p_2, p_1, p_0, p_0ss, p_2s, &faceVolume[face5], faceCentroid[face5], normal[face5]);
28664 computeTri(p_3, p_3s, p_3ss, &faceVolume[face6], faceCentroid[face6], normal[face6]);
28665
28666 computePoly6(p_0ss, p_0s, p_1s, p_3ss, p_3s, p_2s, &area_c, coordinates_c, normalVec_c);
28667
28668 maia::math::vecAvg<3>(M, p_0, p_1, p_2, p_3, p_0s, p_1s, p_2s, p_3s, p_0ss, p_3ss);
28669
28670 for(
MInt i = 0; i < 6; i++) {
28671 computePyra(&faceVolume[*facepointers[i]], faceCentroid[*facepointers[i]], normal[*facepointers[i]], M,
28672 &pyraVolume[i], pyraCentroid[i]);
28673 }
28674 computePyra(&area_c, coordinates_c, normalVec_c, M, &pyraVolume[6], pyraCentroid[6]);
28675 for(
MInt i = 0; i < 7; i++) {
28676 volume_C += pyraVolume[i];
28677 maia::math::vecAdd<3>(coordinates_Cell,
vecScalarMul(pyraVolume[i], pyraCentroid[i], dummy_1));
28678 }
28679 vecScalarMul(F1 / volume_C, coordinates_Cell, coordinates_Cell);
28680
28682
28684 for(
MInt dim = 0; dim < 3; dim++) {
28685 m_bndryCells->
a[bndryId].m_srfcs[0]->m_coordinates[dim] = coordinates_c[dim];
28686 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[dim] = normalVec_c[dim];
28687 m_bndryCells->
a[bndryId].m_coordinates[dim] = coordinates_Cell[dim];
28688 }
28689
28690 break;
28691 }
28692 default: {
28693 stringstream errorMessage;
28694 errorMessage << "FvBndryCnd3D::createCutFaceMGC - Error: Inconsistent type implementation." << endl;
28695 mTerm(1, AT_, errorMessage.str());
28696 }
28697 }
28698
28699 absA = F0;
28700 for(
MInt i = 0; i < nDim; i++) {
28701
28702 if(!nfs_cur[2 * i + 1]) faceVolume[2 * i + 1] = 0;
28703 if(!nfs_cur[2 * i]) faceVolume[2 * i] = 0;
28704
28706
28707 faceDiff[i] = faceVolume[2 * i + 1] - faceVolume[2 * i];
28708 absA +=
POW2(faceDiff[i]);
28709 }
28710
28711 for(
MInt i = 0; i < nDim; i++) {
28712 m_bndryCells->
a[bndryId].m_srfcs[0]->m_normalVector[i] = faceDiff[i] / sqrt(absA);
28713 }
28715
28716
28717
28718 for(
MInt face = 0; face < 6; face++) {
28720 for(
MInt i = 0; i < noFaces; i++) {
28721 if(!(*facepointers[i] == face)) continue;
28722 sideId = face % 2;
28723 spaceId = face / 2;
28724 if(nfs_cur[face]) {
28727 } else {
28731 else
28732 continue;
28733 } else
28734 continue;
28735 }
28737 otherSideId = (sideId + 1) % 2;
28743 for(
MInt dim = 0; dim < nDim; dim++) {
28745 }
28748 break;
28749 } else {
28750 stringstream errorMessage;
28751 errorMessage << " Error in FvBndryCnd3D::createCutFaceMGC - 6. Problem with face " << face << endl;
28752 errorMessage << "bndryId: " << bndryId << endl;
28753 errorMessage <<
"cellId: " <<
cellId << endl;
28754 errorMessage <<
"bndryCnd: " <<
m_bndryCells->
a[bndryId].m_srfcs[0]->m_bndryCndId;
28755 errorMessage << "case, subcase: " << currentCase << ", " << currentSubCase << endl;
28756 errorMessage << "face Volumes: " << faceVolume[0] << ", " << faceVolume[1] << ", " << faceVolume[2] << ", "
28757 << faceVolume[3] << ", " << faceVolume[4] << ", " << faceVolume[5] << endl;
28758 errorMessage << "faceCentroids: ";
28759 for(
MInt j = 0; j < 6; j++) {
28760 errorMessage << faceCentroid[j][0] << ", " << faceCentroid[j][1] << ", " << faceCentroid[j][2] << "; "
28761 << endl;
28762 }
28763 errorMessage << "area_c = " << area_c << endl;
28764 errorMessage << "normalVec_c = " << normalVec_c[0] << ", " << normalVec_c[1] << ", " << normalVec_c[2]
28765 << endl;
28766 errorMessage << "coordinates_c = " << coordinates_c[0] << ", " << coordinates_c[1] << ", "
28767 << coordinates_c[2] << endl;
28768 errorMessage << "coordinates_Cell = " << coordinates_Cell[0] << ", " << coordinates_Cell[1] << ", "
28769 << coordinates_Cell[2] << endl;
28770 errorMessage << "volume_C = " << volume_C << endl;
28771 mTerm(1, AT_, errorMessage.str());
28772 }
28773 }
28774 }
28775
28776 for(
MInt face = 0; face < 6; face++) {
28777 if(!nfs_cur[face]) {
28779 }
28780 }
28781 }
28782 }
28783
28784
28785 if(keepNghbrs) {
28791 }
28797 }
28798 }
28799 }
28800
28801
28802
28803 if(!keepNghbrs) {
28804 for(
MInt ambId = 0; ambId < noAmbiguousCells; ambId++) {
28805 MInt bndryId = ambiguousCells[ambId];
28807
28808
28809 if(splitCell[bndryId] && !splitFace[bndryId]) {
28810
28811 for(
MInt child = 0; child < 3; child++) {
28813 if(bndryId2 == -1) break;
28815
28816
28817 for(
MInt corn = 0; corn < 8; corn++) {
28818 if(cornerCellMapping[ambId * 8 + corn] == cellId2) {
28819 for(
MInt f = 0; f < 3; f++) {
28820 faceTMP = cornerFaceMapping[corn * 3 + f];
28823 if(nghbrBndryId < 0)
mTerm(1, AT_,
"strange split cell neighbor");
28824 if(!splitCell[nghbrBndryId] && !splitFace[nghbrBndryId]) {
28825
28826
28827
28828 srfcId =
m_bndryCells->
a[nghbrBndryId].m_associatedSrfc[opposite[faceTMP]];
28831 } else {
28833 }
28834 } else {
28835 if(splitCell[nghbrBndryId]) {
28836
28837
28838
28839 srfcId =
m_bndryCells->
a[bndryId2].m_associatedSrfc[faceTMP];
28842 } else {
28844 }
28845 } else if(splitFace[nghbrBndryId]) {
28846
28847 }
28848 }
28849 }
28850 } else {
28851
28852 }
28853 }
28854
28855
28856
28857
28858
28859
28860 }
28861
28862
28863 for(
MInt corn = 0; corn < 8; corn++) {
28864 if(cornerCellMapping[ambId * 8 + corn] == cellId) {
28865 for(
MInt f = 0; f < 3; f++) {
28866 faceTMP = cornerFaceMapping[corn * 3 + f];
28869 if(!splitCell[nghbrBndryId] && !splitFace[nghbrBndryId]) {
28870
28871
28872
28873 } else {
28874 if(splitCell[nghbrBndryId]) {
28875
28876
28877
28878 srfcId =
m_bndryCells->
a[bndryId].m_associatedSrfc[faceTMP];
28881 } else {
28883 }
28884 } else if(splitFace[nghbrBndryId]) {
28885
28886 }
28887 }
28888 }
28889 } else {
28890
28891 }
28892 }
28893
28894
28895
28896
28897
28898
28899 }
28900
28901
28902
28903 if(!splitCell[bndryId] && splitFace[bndryId]) {
28904 for(
MInt face = 0; face < 6; face++) {
28906
28907
28908 if(srfcId < -1) {
28909
28912
28915 if(splitCell[nghbrBndryId]) {
28918 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
28919 else
28921 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
28922 } else if(splitFace[nghbrBndryId]) {
28923
28924 } else {
28925
28926 stringstream errorMessage;
28927 errorMessage << " Error in FvBndryCnd3D::createCutFaceMGC::this should not happen! ";
28928 mTerm(1, AT_, errorMessage.str());
28929 }
28930
28931
28934 if(splitCell[nghbrBndryId]) {
28937 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
28938 else
28940 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
28941 } else if(splitFace[nghbrBndryId]) {
28942
28943 } else {
28944
28945 }
28946 } else if(srfcId >= 0) {
28947
28950 if(splitCell[nghbrBndryId]) {
28951 for(
MInt i = 0; i < 4; i++) {
28952 if(cornerCellMapping[ambIds[nghbrBndryId] * 8 + faceCornerMapping[face * 6 + i]] == -1) {
28953
28954 } else {
28955
28956
28957
28960 else
28962 break;
28963 }
28964 }
28965 } else if(splitFace[nghbrBndryId]) {
28966
28967 } else {
28968
28969 }
28970 }
28971 }
28972 }
28973 }
28974
28975
28976 } else {
28977 for(
MInt ambId = 0; ambId < noAmbiguousCells; ambId++) {
28978 MInt bndryId = ambiguousCells[ambId];
28980
28981
28982 if(splitCell[bndryId] && !splitFace[bndryId]) {
28983
28984 for(
MInt child = 0; child < 3; child++) {
28986 if(bndryId2 == -1) break;
28988 MInt splitCellId = cellId2;
28989 for(
MInt corn = 0; corn < 8; corn++) {
28990 if(cornerCellMapping[ambId * 8 + corn] == cellId2) {
28991 for(
MInt f = 0; f < 3; f++) {
28992 faceTMP = cornerFaceMapping[corn * 3 + f];
28995 }
28999 if(!splitCell[nghbrBndryId] && !splitFace[nghbrBndryId]) {
29000
29001
29002
29005 srfcId =
m_bndryCells->
a[nghbrBndryId].m_associatedSrfc[opposite[faceTMP]];
29006 if(srfcId > -1) {
29009 } else {
29011 }
29012 } else
29015 } else {
29016 if(splitCell[nghbrBndryId]) {
29017
29019 cornerCellMapping[ambIds[nghbrBndryId] * 8 + corn + neighborCorner[faceTMP]];
29020 srfcId =
m_bndryCells->
a[bndryId2].m_associatedSrfc[faceTMP];
29021 if(srfcId > -1) {
29025 } else {
29028 }
29029 }
29030 } else if(splitFace[nghbrBndryId]) {
29031
29032
29034 }
29035 }
29036 }
29037 } else {
29038
29039 }
29040 }
29041
29044 }
29045 }
29046
29047
29048 for(
MInt corn = 0; corn < 8; corn++) {
29049 if(cornerCellMapping[ambId * 8 + corn] == cellId) {
29050 for(
MInt f = 0; f < 3; f++) {
29051 faceTMP = cornerFaceMapping[corn * 3 + f];
29054 if(!splitCell[nghbrBndryId] && !splitFace[nghbrBndryId]) {
29055
29056
29057
29058
29061 } else {
29062 if(splitCell[nghbrBndryId]) {
29063
29065 cornerCellMapping[ambIds[nghbrBndryId] * 8 + corn + neighborCorner[faceTMP]];
29066 srfcId =
m_bndryCells->
a[bndryId].m_associatedSrfc[faceTMP];
29069 } else {
29071 }
29072 } else if(splitFace[nghbrBndryId]) {
29074
29075 }
29076 }
29077 }
29078 } else {
29079
29080 }
29081 }
29082
29085 }
29086 }
29087
29088
29089
29090 if(!splitCell[bndryId] && splitFace[bndryId]) {
29091 for(
MInt face = 0; face < 6; face++) {
29093
29094
29095 if(srfcId < -1) {
29096
29099
29102 if(splitCell[nghbrBndryId]) {
29104 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
29107 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
29108 else
29110 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
29111 } else if(splitFace[nghbrBndryId]) {
29112
29114 } else {
29115
29116 stringstream errorMessage;
29117 errorMessage << " Error in FvBndryCnd3D::createCutFaceMGC::this should not happen! ";
29118 mTerm(1, AT_, errorMessage.str());
29119 }
29120
29121
29124 if(splitCell[nghbrBndryId]) {
29125
29127 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
29130 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
29131 else
29133 cornerCellMapping[ambIds[nghbrBndryId] * 8 + splitCorner2];
29134 } else if(splitFace[nghbrBndryId]) {
29135
29137 } else {
29138
29139 stringstream errorMessage;
29140 errorMessage << " Error in FvBndryCnd3D::createCutFaceMGC::this should not happen! ";
29141 mTerm(1, AT_, errorMessage.str());
29142 }
29143 } else if(srfcId >= 0) {
29144
29147 if(splitCell[nghbrBndryId]) {
29148 for(
MInt i = 0; i < 4; i++) {
29149 if(cornerCellMapping[ambIds[nghbrBndryId] * 8 + faceCornerMapping[face * 6 + i]] == -1) {
29150
29151 } else {
29152
29154 cornerCellMapping[ambIds[nghbrBndryId] * 8 + faceCornerMapping[face * 6 + i]];
29157 else
29159 break;
29160 }
29161 }
29162 } else {
29163
29165 }
29166 }
29167 }
29168 }
29169 }
29170
29171
29172 if(keepNghbrs) {
29179 }
29180 }
29181 }
29182 }
29183 }
29184
29185#ifndef NDEBUG
29186 for(
MInt i = 0; i < 15; i++) {
29187 cerr << "Occurences case " << i << " : " << presentCases[i] << " times" << endl;
29188 }
29189 cerr << "---------------------------------------------------------------" << endl << endl << endl;
29191#endif
29192}
void computeTetra(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
MFloat * vecScalarMul(MFloat, MFloat *, MFloat *)
void computeTri(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
void computePoly5(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
void correctNormal(MFloat *)
void correctCell(MFloat *, MFloat *, MFloat *, MFloat *)
void correctFace(MFloat *, MFloat *, MFloat *, MFloat *)
void computePoly4(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
void computePoly6(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
MInt createSplitCell_MGC(MInt, MInt)
produces an exact copy of a fvcell and the respective boundary cell
void computeTrapez(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
void computePyra(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
void computePoly3(MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
void assertValidGridCellId(const MInt cellId) const
Cecks wether the cell cellId is an actual grid-cell.
T * getPointer() const
Deprecated: use begin() instead!
pointer p
Deprecated: use [] instead!
1.9.5