diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp index 45e19235bc3dbe96137342d337043c014a53ba4a..fd796008a73553cffef386033c8f4db2f6554ba2 100644 --- a/Common/Gmsh.cpp +++ b/Common/Gmsh.cpp @@ -244,7 +244,7 @@ int GmshBatch() solver_batch_cb(0, (void*)CTX::instance()->launchSolverAtStartup); #endif - if (false) { + if (true) { // 11/06/13 : This will be removed later ! for (int i = 1; i < MSH_NUM_TYPE+1; ++i) { if (i == 76 || i == 77 || i == 78) @@ -253,24 +253,24 @@ int GmshBatch() Msg::Warning("ignoring unknown %d", i); continue; } - if (MElement::ParentTypeFromTag(i) == TYPE_PRI) { - Msg::Info("ignoring prism %d", i); + if (MElement::ParentTypeFromTag(i) == TYPE_PRI && MElement::OrderFromTag(i) > 2) { + Msg::Warning("ignoring prism %d (different implementation)", i); continue; } if (MElement::ParentTypeFromTag(i) == TYPE_PYR) { - Msg::Info("ignoring pyramid %d", i); + Msg::Warning("ignoring pyramid %d", i); continue; } if (MElement::ParentTypeFromTag(i) == TYPE_POLYG) { - Msg::Info("ignoring polygone %d", i); + Msg::Warning("ignoring polygone %d", i); continue; } if (MElement::ParentTypeFromTag(i) == TYPE_POLYH) { - Msg::Info("ignoring polyh�dre %d", i); + Msg::Warning("ignoring polyh�dre %d", i); continue; } if (MElement::ParentTypeFromTag(i) == TYPE_XFEM) { - Msg::Info("ignoring xfem %d", i); + Msg::Warning("ignoring xfem %d", i); continue; } const nodalBasis *fs = BasisFactory::getNodalBasis(i); diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp index 5f12e1aac5997018ee59653e5cddb605d704d39d..394720ec305f16ed6693d8e3e607e209c52d1def 100644 --- a/Geo/MElement.cpp +++ b/Geo/MElement.cpp @@ -1262,6 +1262,20 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name) case MSH_PRI_6 : if(name) *name = "Prism 6"; return 6; case MSH_PRI_15 : if(name) *name = "Prism 15"; return 6 + 9; case MSH_PRI_18 : if(name) *name = "Prism 18"; return 6 + 9 + 3; + case MSH_PRI_40 : if(name) *name = "Prism 40"; return 6 + 18 + 12+2 + 2*1; + case MSH_PRI_75 : if(name) *name = "Prism 75"; return 6 + 27 + 27+6 + 3*3; + case MSH_PRI_126 : if(name) *name = "Prism 126"; return 6 + 36 + 48+12 + 4*6; + case MSH_PRI_196 : if(name) *name = "Prism 196"; return 6 + 45 + 75+20 + 5*10; + case MSH_PRI_288 : if(name) *name = "Prism 288"; return 6 + 54 + 108+30 + 6*15; + case MSH_PRI_405 : if(name) *name = "Prism 405"; return 6 + 63 + 147+42 + 7*21; + case MSH_PRI_550 : if(name) *name = "Prism 550"; return 6 + 72 + 196+56 + 8*28; + case MSH_PRI_38 : if(name) *name = "Prism 38"; return 6 + 18 + 12+2; + case MSH_PRI_66 : if(name) *name = "Prism 66"; return 6 + 27 + 27+6; + case MSH_PRI_102 : if(name) *name = "Prism 102"; return 6 + 36 + 48+12; + case MSH_PRI_146 : if(name) *name = "Prism 146"; return 6 + 45 + 75+20; + case MSH_PRI_198 : if(name) *name = "Prism 198"; return 6 + 54 + 108+30; + case MSH_PRI_258 : if(name) *name = "Prism 258"; return 6 + 63 + 147+42; + case MSH_PRI_326 : if(name) *name = "Prism 326"; return 6 + 72 + 196+56; case MSH_PYR_1 : if(name) *name = "Pyramid 1"; return 1; case MSH_PYR_5 : if(name) *name = "Pyramid 5"; return 5; case MSH_PYR_13 : if(name) *name = "Pyramid 13"; return 5 + 8; diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp index 36730f595046f308653cfaa78cc2de83e15f9a4c..9f15a203d45072af512d4df18c39292b7114c82d 100644 --- a/Numeric/BasisFactory.cpp +++ b/Numeric/BasisFactory.cpp @@ -40,7 +40,7 @@ const nodalBasis* BasisFactory::getNodalBasis(int elementType) F = new pyramidalBasis(elementType); break; default: - Msg::Error("Unknown type of element."); + Msg::Error("Unknown type of element %d (in BasisFactory)", elementType); return NULL; } diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp index bd3411536c52b01b90accf83090b0753a7cd6a5f..21a69e2eb5dbcf7960ecbc44dfe8de5907dc7022 100644 --- a/Numeric/nodalBasis.cpp +++ b/Numeric/nodalBasis.cpp @@ -20,6 +20,7 @@ int nodalBasis::compareNewAlgoPointsWithOld() const } if (points_newAlgo.size1() != points.size1()) { Msg::Error("%d: pas meme taille (%d, %d, %d) %s", type, parentType, order, serendip, *name); + Msg::Error(" |--> points[%d, %d] vs newPoints[%d, %d]", points.size1(), points.size2(), points_newAlgo.size1(), points_newAlgo.size2()); return 2; } for (int i = 0; i < points.size1(); ++i) { @@ -1247,8 +1248,8 @@ static void generateFaceClosureHexFull(nodalBasis::clCont &closure, static void getFaceClosurePrism(int iFace, int iSign, int iRotate, nodalBasis::closure &closure, int order) { - if (order > 2) - Msg::Error("FaceClosure not implemented for prisms of order %d",order); + //if (order > 2) + // Msg::Error("FaceClosure not implemented for prisms of order %d",order); bool isTriangle = iFace<2; int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1); closure.clear(); @@ -1285,6 +1286,8 @@ static void getFaceClosurePrism(int iFace, int iSign, int iRotate, static void generateFaceClosurePrism(nodalBasis::clCont &closure, int order) { + if (order > 2) + Msg::Error("FaceClosure not implemented for prisms of order %d",order); closure.clear(); for (int iRotate = 0; iRotate < 4; iRotate++){ for (int iSign = 1; iSign >= -1; iSign -= 2){ @@ -1603,7 +1606,7 @@ nodalBasis::nodalBasis(int tag) closureRef.resize(48, 0); } else { - generateFaceClosurePrism(closures, order); + if (order < 2) generateFaceClosurePrism(closures, order); generateFaceClosurePrismFull(fullClosures, closureRef, order); } break; diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp index bb50681b9197417236d51cbbdb574d100765aad3..5ae836323f1400b65c0b7fd3808dd3d6ead67b16 100644 --- a/Numeric/pointsGenerators.cpp +++ b/Numeric/pointsGenerators.cpp @@ -8,6 +8,7 @@ #include "MQuadrangle.h" #include "MTetrahedron.h" #include "MHexahedron.h" +#include "MPrism.h" /* --- Lines --- */ @@ -832,9 +833,9 @@ fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip) return monomials; } -fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order) +fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints) { - int nbMonomials = (order+1)*(order+1); + int nbMonomials = forSerendipPoints ? order*4 : (order+1)*(order+1); fullMatrix<int> monomials(nbMonomials, 2); monomials(0, 0) = 0; @@ -864,9 +865,12 @@ fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order) monomials(index, 1) = monomials(i0, 1) + u_1 * i; } } - fullMatrix<int> inner = gmshGenerateMonomialsQuadrangle(order-2); - inner.add(1); - monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0); + + if (!forSerendipPoints) { + fullMatrix<int> inner = gmshGenerateMonomialsQuadrangle(order-2); + inner.add(1); + monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0); + } } } return monomials; @@ -956,59 +960,115 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip) return monomials; } -/* -fullMatrix<int> gmshGenerateMonomialsPrism(int order) +fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints) { - const int prism18Pts[18][3] = { - {0, 0, 0}, // 0 - {2, 0, 0}, // 1 - {0, 2, 0}, // 2 - {0, 0, 2}, // 3 - {2, 0, 2}, // 4 - {0, 2, 2}, // 5 - {1, 0, 0}, // 6 - {0, 1, 0}, // 7 - {0, 0, 1}, // 8 - {1, 1, 0}, // 9 - {2, 0, 1}, // 10 - {0, 2, 1}, // 11 - {1, 0, 2}, // 12 - {0, 1, 2}, // 13 - {1, 1, 2}, // 14 - {1, 0, 1}, // 15 - {0, 1, 1}, // 16 - {1, 1, 1}, // 17 - }; - - int nbMonomials = (order + 1)*(order + 1)*(order + 2)/2; + int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 : (order + 1) * (order + 1)*(order + 2)/2; fullMatrix<int> monomials(nbMonomials, 3); - int index = 0; - fullMatrix<int> triMonomials = gmshGenerateMonomialsTriangle(order,false); - fullMatrix<int> lineMonomials = gmshGenerateMonomialsLine(order); + monomials(0, 0) = 0; + monomials(0, 1) = 0; + monomials(0, 2) = 0; - if (order == 2) - for (int i =0; i<18; i++) - for (int j=0; j<3;j++) - monomials(i,j) = prism18Pts[i][j]; - else - for (int j = 0; j <lineMonomials.size1() ; j++) { - for (int i = 0; i < triMonomials.size1(); i++) { - monomials(index,0) = triMonomials(i,0); - monomials(index,1) = triMonomials(i,1); - monomials(index,2) = lineMonomials(j,0); - index++; + if (order > 0) { + monomials(1, 0) = order; + monomials(1, 1) = 0; + monomials(1, 2) = 0; + + monomials(2, 0) = 0; + monomials(2, 1) = order; + monomials(2, 2) = 0; + + monomials(3, 0) = 0; + monomials(3, 1) = 0; + monomials(3, 2) = order; + + monomials(4, 0) = order; + monomials(4, 1) = 0; + monomials(4, 2) = order; + + monomials(5, 0) = 0; + monomials(5, 1) = order; + monomials(5, 2) = order; + + if (order > 1) { + int index = 6; + for (int iedge = 0; iedge < 9; ++iedge) { + int i0 = MPrism::edges_prism(iedge, 0); + int i1 = MPrism::edges_prism(iedge, 1); + + int u_1 = (monomials(i1,0)-monomials(i0,0)) / order; + int u_2 = (monomials(i1,1)-monomials(i0,1)) / order; + int u_3 = (monomials(i1,2)-monomials(i0,2)) / order; + + for (int i = 1; i < order; ++i, ++index) { + monomials(index, 0) = monomials(i0, 0) + i * u_1; + monomials(index, 1) = monomials(i0, 1) + i * u_2; + monomials(index, 2) = monomials(i0, 2) + i * u_3; + } } - } - return monomials; + if (!forSerendipPoints) { + fullMatrix<int> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2); + dudvQ.add(1); + + fullMatrix<int> dudvT; + if (order > 2) dudvT = gmshGenerateMonomialsTriangle(order - 3); + dudvT.add(1); + + for (int iface = 0; iface < 5; ++iface) { + int i0, i1, i2; + i0 = MPrism::faces_prism(iface, 0); + i1 = MPrism::faces_prism(iface, 1); + fullMatrix<int> dudv; + if (MPrism::faces_prism(iface, 3) != -1) { + i2 = MPrism::faces_prism(iface, 3); + dudv.setAsProxy(dudvQ); + } + else if (order > 2) { + i2 = MPrism::faces_prism(iface, 2); + dudv.setAsProxy(dudvT); + } + else { + continue; + } + + int u[3]; + u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order; + u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order; + u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order; + int v[3]; + v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order; + v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order; + v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order; + + for (int i = 0; i < dudv.size1(); ++i, ++index) { + monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1); + monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1); + monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1); + } + } + + if (order > 2) { + fullMatrix<int> triMonomials = gmshGenerateMonomialsTriangle(order - 3); + fullMatrix<int> lineMonomials = gmshGenerateMonomialsLine(order - 2); + for (int i = 0; i < triMonomials.size1(); ++i) { + for (int j = 0; j < lineMonomials.size1(); ++j, ++index) { + monomials(index, 0) = 1 + triMonomials(i, 0); + monomials(index, 1) = 1 + triMonomials(i, 1); + monomials(index, 2) = 1 + lineMonomials(j, 0); + } + } + } + } + } + } + return monomials; } -*/ -fullMatrix<int> gmshGenerateMonomialsHexahedron(int order) +fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints) { - int nbMonomials = (order+1)*(order+1)*(order+1); + int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 : (order+1)*(order+1)*(order+1); fullMatrix<int> monomials(nbMonomials, 3); monomials(0, 0) = 0; @@ -1064,30 +1124,32 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order) fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 2); dudv.add(1); - for (int iface = 0; iface < 6; ++iface) { - int i0 = MHexahedron::faces_hexa(iface, 0); - int i1 = MHexahedron::faces_hexa(iface, 1); - int i3 = MHexahedron::faces_hexa(iface, 3); + if (!forSerendipPoints) { + for (int iface = 0; iface < 6; ++iface) { + int i0 = MHexahedron::faces_hexa(iface, 0); + int i1 = MHexahedron::faces_hexa(iface, 1); + int i3 = MHexahedron::faces_hexa(iface, 3); - int u[3]; - u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order; - u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order; - u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order; - int v[3]; - v[0] = (monomials(i3, 0) - monomials(i0, 0)) / order; - v[1] = (monomials(i3, 1) - monomials(i0, 1)) / order; - v[2] = (monomials(i3, 2) - monomials(i0, 2)) / order; - - for (int i = 0; i < dudv.size1(); ++i, ++index) { - monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1); - monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1); - monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1); + int u[3]; + u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order; + u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order; + u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order; + int v[3]; + v[0] = (monomials(i3, 0) - monomials(i0, 0)) / order; + v[1] = (monomials(i3, 1) - monomials(i0, 1)) / order; + v[2] = (monomials(i3, 2) - monomials(i0, 2)) / order; + + for (int i = 0; i < dudv.size1(); ++i, ++index) { + monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1); + monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1); + monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1); + } } - } - fullMatrix<int> inner = gmshGenerateMonomialsHexahedron(order - 2); - inner.add(1); - monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0); + fullMatrix<int> inner = gmshGenerateMonomialsHexahedron(order - 2); + inner.add(1); + monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0); + } } } return monomials; diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h index 8199e874f03dc61ad199070e61f180c16adf740e..4524266ac9288c2cf14d727818128dd20bd86397 100644 --- a/Numeric/pointsGenerators.h +++ b/Numeric/pointsGenerators.h @@ -35,12 +35,14 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip); fullMatrix<int> gmshGenerateMonomialsLine(int order); fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip = false); -fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order); -//fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order); +fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints = false); +//fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order); //FIXME fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip = false); -fullMatrix<int> gmshGenerateMonomialsHexahedron(int order); -//fullMatrix<int> gmshGenerateMonomialsPrism(int order); +fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints = false); +//fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order); //FIXME +fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false); +//fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order); //FIXME //fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool serendip); diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp index bdecd24afdbb8a796540cdf03701ce09f92aa566..d4cef22282809c3f54f24709ec9ba66c3f8f83af 100644 --- a/Numeric/polynomialBasis.cpp +++ b/Numeric/polynomialBasis.cpp @@ -388,40 +388,57 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag) // TEST NEW ALGO POINTS / MONOMIALS - bool centered = false; + int rescale = 0; switch (parentType) { case TYPE_PNT : monomials_newAlgo = gmshGenerateMonomialsLine(0); break; case TYPE_LIN : monomials_newAlgo = gmshGenerateMonomialsLine(order); - centered = true; + rescale = 1; break; case TYPE_TRI : monomials_newAlgo = gmshGenerateMonomialsTriangle(order, serendip); break; case TYPE_QUA : - if (serendip) return; - monomials_newAlgo = gmshGenerateMonomialsQuadrangle(order); - centered = true; + monomials_newAlgo = gmshGenerateMonomialsQuadrangle(order, serendip); + rescale = 1; break; case TYPE_TET : monomials_newAlgo = gmshGenerateMonomialsTetrahedron(order, serendip); break; case TYPE_HEX : - if (serendip) return; - monomials_newAlgo = gmshGenerateMonomialsHexahedron(order); - centered = true; + monomials_newAlgo = gmshGenerateMonomialsHexahedron(order, serendip); + rescale = 1; + break; + case TYPE_PRI : + monomials_newAlgo = gmshGenerateMonomialsPrism(order, serendip); + rescale = 2; break; } copy(monomials_newAlgo, points_newAlgo); if (order == 0) return; - if (centered) { - points_newAlgo.scale(2./order); - points_newAlgo.add(-1.); - } - else { - points_newAlgo.scale(1./order); + switch (rescale) { + case 0 : + points_newAlgo.scale(1./order); + break; + case 1 : + points_newAlgo.scale(2./order); + points_newAlgo.add(-1.); + break; + case 2 : + { + fullMatrix<double> tmp; + tmp.setAsProxy(points_newAlgo, 0, 2); + tmp.scale(1./order); + + tmp.setAsProxy(points_newAlgo, 2, 1); + tmp.scale(2./order); + tmp.add(-1.); + break; + } + default : + break; } }