diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index bc79a563..40a399eb 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -43,8 +43,8 @@ namespace netgen bcnames.SetSize(0); cd2names.SetSize(0); + this->comm = netgen :: ng_comm; #ifdef PARALLEL - this->comm = MPI_COMM_WORLD; paralleltop = new ParallelMeshTopology (*this); #endif } @@ -83,6 +83,10 @@ namespace netgen #endif } + void Mesh :: SetCommunicator(MPI_Comm acomm) + { + this->comm = acomm; + } Mesh & Mesh :: operator= (const Mesh & mesh2) { @@ -314,7 +318,7 @@ namespace netgen SurfaceElementIndex si = surfelements.Size(); surfelements.Append (el); - if (el.index > facedecoding.Size()) + if (el.index<=0 || el.index > facedecoding.Size()) cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl; surfelements.Last().next = facedecoding[el.index-1].firstelement; @@ -1173,7 +1177,7 @@ namespace netgen for (i = 1; i<= GetNSeg(); i++) { Segment & seg = LineSegment(i); - if ( seg.cd2i <= n ) + if ( seg.edgenr <= n ) seg.SetBCName (cd2names[seg.edgenr-1]); else seg.SetBCName(0); @@ -1304,7 +1308,7 @@ namespace netgen } - void Mesh :: DoArchive (ngstd::Archive & archive) + void Mesh :: DoArchive (Archive & archive) { archive & dimension; archive & points; @@ -1316,25 +1320,14 @@ namespace netgen archive & *ident; - - // archive geometry - if (archive.Output()) - { - ostringstream ost; - if (geometry) - geometry -> SaveToMeshFile (ost); - archive << ost.str(); - } - else - { - string str; - archive & str; - istringstream ist(str); - geometry = geometryregister.LoadFromMeshFile (ist); - } + archive.Shallow(geometry); + archive & *curvedelems; if (archive.Input()) { + int rank = MyMPI_GetId(GetCommunicator()); + int ntasks = MyMPI_GetNTasks(GetCommunicator()); + RebuildSurfaceElementLists(); CalcSurfacesOfNode (); @@ -1687,6 +1680,7 @@ namespace netgen void Mesh :: CalcSurfacesOfNode () { + static Timer t("Mesh::CalcSurfacesOfNode"); RegionTimer reg (t); // surfacesonnode.SetSize (GetNP()); TABLE surfacesonnode(GetNP()); @@ -1871,8 +1865,7 @@ namespace netgen void Mesh :: FindOpenElements (int dom) { - static int timer = NgProfiler::CreateTimer ("Mesh::FindOpenElements"); - NgProfiler::RegionTimer reg (timer); + static Timer t("Mesh::FindOpenElements"); RegionTimer reg (t); int np = GetNP(); int ne = GetNE(); @@ -2724,6 +2717,8 @@ namespace netgen void Mesh :: CalcLocalH (double grading) { + static Timer t("Mesh::CalcLocalH"); RegionTimer reg(t); + if (!lochfunc) { Point3d pmin, pmax; @@ -3249,6 +3244,8 @@ namespace netgen void Mesh :: Compress () { + static Timer t("Mesh::Compress"); RegionTimer reg(t); + Array op2np(GetNP()); Array hpoints; BitArrayChar pused(GetNP()); @@ -3388,7 +3385,7 @@ namespace netgen for (int i = 0; i < lockedpoints.Size(); i++) lockedpoints[i] = op2np[lockedpoints[i]]; - + /* for (int i = 0; i < facedecoding.Size(); i++) facedecoding[i].firstelement = -1; for (int i = surfelements.Size()-1; i >= 0; i--) @@ -3397,7 +3394,8 @@ namespace netgen surfelements[i].next = facedecoding[ind-1].firstelement; facedecoding[ind-1].firstelement = i; } - + */ + RebuildSurfaceElementLists (); CalcSurfacesOfNode(); @@ -3515,6 +3513,8 @@ namespace netgen int Mesh :: CheckOverlappingBoundary () { + static Timer t("Mesh::CheckOverlappingBoundary"); RegionTimer reg(t); + int i, j, k; Point3d pmin, pmax; @@ -4327,10 +4327,12 @@ namespace netgen elementsearchtree = NULL; int ne = (dimension == 2) ? GetNSE() : GetNE(); - + if (dimension == 3 && !GetNE() && GetNSE()) + ne = GetNSE(); + if (ne) { - if (dimension == 2) + if (dimension == 2 || (dimension == 3 && !GetNE()) ) { Box<3> box (Box<3>::EMPTY_BOX); for (SurfaceElementIndex sei = 0; sei < ne; sei++) @@ -4393,7 +4395,11 @@ namespace netgen sol.Y() = (-a12 * aTrhs.X() + a11 * aTrhs.Y()) / det; return 0; } - + + bool ValidBarCoord(double lami[3], double eps=1e-12) + { + return (lami[0]<=1.+eps && lami[0]>=0.-eps && lami[1]<=1.+eps && lami[1]>=0.-eps && lami[2]<=1.+eps && lami[2]>=0.-eps ); + } bool Mesh :: PointContainedIn2DElement(const Point3d & p, double lami[3], @@ -4424,16 +4430,167 @@ namespace netgen Vec3d c = p4 - a; Vec3d d = p3 - a - b - c; - + /*cout << "p = " << p << endl; + cout << "p1 = " << p1 << endl; + cout << "p2 = " << p2 << endl; + cout << "p3 = " << p3 << endl; + cout << "p4 = " << p4 << endl; + + cout << "a = " << a << endl; + cout << "b = " << b << endl; + cout << "c = " << c << endl; + cout << "d = " << d << endl;*/ + + + Vec3d pa = p-a; double dxb = d.X()*b.Y()-d.Y()*b.X(); - double dxc = d.X()*c.Y()-d.Y()*c.X(); + double dxc = d.X()*c.Y()-d.Y()*c.X(); + double bxc = b.X()*c.Y()-b.Y()*c.X(); + double bxpa = b.X()*pa.Y()-b.Y()*pa.X(); + double cxpa = c.X()*pa.Y()-c.Y()*pa.X(); + double dxpa = d.X()*pa.Y()-d.Y()*pa.X(); + + /*cout << "dxb = " << dxb << endl; + cout << "dxc = " << dxc << endl; + cout << "bxc = " << bxc << endl; + cout << "bxpa = " << bxpa << endl; + cout << "cxpa = " << cxpa << endl; + cout << "dxpa = " << dxpa << endl;*/ + + /* + P = a + b x + c y + d x y + 1) P1 = a1 + b1 x + c1 y + d1 x y + 2) P2 = a2 + b2 x + c2 y + d2 x y + + -> det(x,d) = det(a,d) + det(b,d) x + det(c,d) y + -> x = 1/det(b,d) *( det(P-a,d)-det(c,d) y ) + -> y = 1/det(c,d) *( det(P-a,d)-det(b,d) x ) + + -> x = (P1 - a1 - c1 y)/(b1 + d1 y) + -> det(c,d) y**2 + [det(d,P-a) + det(c,b)] y + det(b,P-a) = 0 + ( same if we express x = (P2 - a2 - c2 y)/(b2 + d2 y) ) + + -> y = (P1 - a1 - b1 x)/(c1 + d1 x) + -> det(b,d) x**2 + [det(d,P-a) + det(b,c)] x + det(c,P-a) = 0 + ( same if we express y = (P2 - a2 - b2 x)/(c2 + d2 x) + */ + + lami[2]=0.; + double eps = 1.E-12; + double c1,c2,r; + + //First check if point is "exactly" a vertex point + Vec3d d1 = p-p1; + Vec3d d2 = p-p2; + Vec3d d3 = p-p3; + Vec3d d4 = p-p4; + + //cout << " d1 = " << d1 << ", d2 = " << d2 << ", d3 = " << d3 << ", d4 = " << d4 << endl; + + if (d1.Length2() < sqr(eps)*d2.Length2() && d1.Length2() < sqr(eps)*d3.Length2() && d1.Length2() < sqr(eps)*d4.Length2()) + { + lami[0] = lami[1] = 0.; + return true; + } + else if (d2.Length2() < sqr(eps)*d1.Length2() && d2.Length2() < sqr(eps)*d3.Length2() && d2.Length2() < sqr(eps)*d4.Length2()) + { + lami[0] = 1.; + lami[1] = 0.; + return true; + } + else if (d3.Length2() < sqr(eps)*d1.Length2() && d3.Length2() < sqr(eps)*d2.Length2() && d3.Length2() < sqr(eps)*d4.Length2()) + { + lami[0] = lami[1] = 1.; + return true; + } + else if (d4.Length2() < sqr(eps)*d1.Length2() && d4.Length2() < sqr(eps)*d2.Length2() && d4.Length2() < sqr(eps)*d3.Length2()) + { + lami[0] = 0.; + lami[1] = 1.; + return true; + }//if d is nearly 0: solve resulting linear system + else if (d.Length2() < sqr(eps)*b.Length2() && d.Length2() < sqr(eps)*c.Length2()) + { + Vec2d sol; + SolveLinearSystemLS (b, c, p-a, sol); + lami[0] = sol.X(); + lami[1] = sol.Y(); + return ValidBarCoord(lami, eps); + }// if dxc is nearly 0: solve resulting linear equation for y and compute x + else if (fabs(dxc) < sqr(eps)) + { + lami[1] = -bxpa/(dxpa-bxc); + lami[0] = (dxpa-dxc*lami[1])/dxb; + return ValidBarCoord(lami, eps); + }// if dxb is nearly 0: solve resulting linear equation for x and compute y + else if (fabs(dxb) < sqr(eps)) + { + lami[0] = -cxpa/(dxpa+bxc); + lami[1] = (dxpa-dxb*lami[0])/dxc; + return ValidBarCoord(lami, eps); + }//if dxb >= dxc: solve quadratic equation in y and compute x + else if (fabs(dxb) >= fabs(dxc)) + { + c1 = (bxc-dxpa)/dxc; + c2 = -bxpa/dxc; + r = c1*c1/4.0-c2; + + //quadratic equation has only 1 (unstable) solution + if (fabs(r) < eps) //not eps^2! + { + lami[1] = -c1/2; + lami[0] = (dxpa-dxc*lami[1])/dxb; + return ValidBarCoord(lami, eps); + } + if (r < 0) return false; + + lami[1] = -c1/2+sqrt(r); + lami[0] = (dxpa-dxc*lami[1])/dxb; + + if (ValidBarCoord(lami, eps)) + return true; + else + { + lami[1] = -c1/2-sqrt(r); + lami[0] = (dxpa-dxc*lami[1])/dxb; + return ValidBarCoord(lami, eps); + } + }//if dxc > dxb: solve quadratic equation in x and compute y + else + { + c1 = (-bxc-dxpa)/dxb; + c2 = -cxpa/dxb; + r = c1*c1/4.0-c2; + + //quadratic equation has only 1 (unstable) solution + if (fabs(r) < eps) //not eps^2! + { + lami[0] = -c1/2; + lami[1] = (dxpa-dxb*lami[0])/dxc; + return ValidBarCoord(lami, eps); + } + if (r < 0) return false; + + lami[0] = -c1/2+sqrt(r); + lami[1] = (dxpa-dxb*lami[0])/dxc; + + if (ValidBarCoord(lami, eps)) + return true; + else + { + lami[0] = -c1/2-sqrt(r); + lami[1] = (dxpa-dxb*lami[0])/dxc; + return ValidBarCoord(lami, eps); + } + } + + /* double dxa = d.X()*a.Y()-d.Y()*a.X(); double dxp = d.X()*p.Y()-d.Y()*p.X(); double c0,c1,c2; // ,rt; - lami[2]=0.; - double eps = 1.E-12; + Vec3d dp13 = p3-p1; Vec3d dp24 = p4-p2; @@ -4453,12 +4610,12 @@ namespace netgen if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) return true; - /* - lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/ - (b.X()*c.Y() -b.Y()*c.X()); - lami[1]=(-b.Y()*(p.X()-a.X())+b.X()*(p.Y()-a.Y()))/ - (b.X()*c.Y() -b.Y()*c.X()); - */ + + //lami[0]=(c.Y()*(p.X()-a.X())-c.X()*(p.Y()-a.Y()))/ + //(b.X()*c.Y() -b.Y()*c.X()); + //lami[1]=(-b.Y()*(p.X()-a.X())+b.X()*(p.Y()-a.Y()))/ + // (b.X()*c.Y() -b.Y()*c.X()); + } else if(fabs(dxb) <= eps*fabs(dxc)) @@ -4552,7 +4709,7 @@ namespace netgen if(lami[1]<=1.+eps && lami[1]>=0.-eps && lami[0]<=1.+eps && lami[0]>=0.-eps) return true; - } + }*/ //cout << "lam0,1 = " << lami[0] << ", " << lami[1] << endl; @@ -4712,7 +4869,7 @@ namespace netgen lam(2) > -eps && lam(0) + lam(1) + lam(2) < 1+eps); } - else if (el.GetType() == PRISM) + else if (el.GetType() == PRISM || el.GetType() == PRISM15) { retval = (lam(0) > -eps && lam(1) > -eps && @@ -4720,7 +4877,7 @@ namespace netgen lam(2) < 1+eps && lam(0) + lam(1) < 1+eps); } - else if (el.GetType() == PYRAMID) + else if (el.GetType() == PYRAMID || el.GetType() == PYRAMID13) { retval = (lam(0) > -eps && lam(1) > -eps && @@ -4728,7 +4885,7 @@ namespace netgen lam(0) + lam(2) < 1+eps && lam(1) + lam(2) < 1+eps); } - else if (el.GetType() == HEX) + else if (el.GetType() == HEX || el.GetType() == HEX20) { retval = (lam(0) > -eps && lam(0) < 1+eps && lam(1) > -eps && lam(1) < 1+eps && @@ -4831,10 +4988,7 @@ namespace netgen bool build_searchtree, const bool allowindex) const { - const double pointtol = 1e-12; - netgen::Point<3> pmin = p - Vec<3> (pointtol, pointtol, pointtol); - netgen::Point<3> pmax = p + Vec<3> (pointtol, pointtol, pointtol); - if (dimension == 2) + if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE())) { int ne; int ps_startelement = 0; // disable global buffering @@ -4846,8 +5000,11 @@ namespace netgen if (elementsearchtree || build_searchtree) { // update if necessary: - const_cast(*this).BuildElementSearchTree (); - elementsearchtree->GetIntersecting (pmin, pmax, locels); + const_cast(*this).BuildElementSearchTree (); + // double tol = elementsearchtree->Tolerance(); + // netgen::Point<3> pmin = p - Vec<3> (tol, tol, tol); + // netgen::Point<3> pmax = p + Vec<3> (tol, tol, tol); + elementsearchtree->GetIntersecting (p, p, locels); ne = locels.Size(); } else @@ -4889,8 +5046,11 @@ namespace netgen if (elementsearchtree || build_searchtree) { // update if necessary: - const_cast(*this).BuildElementSearchTree (); - elementsearchtree->GetIntersecting (pmin, pmax, locels); + const_cast(*this).BuildElementSearchTree (); + // double tol = elementsearchtree->Tolerance(); + // netgen::Point<3> pmin = p - Vec<3> (tol, tol, tol); + // netgen::Point<3> pmax = p + Vec<3> (tol, tol, tol); + elementsearchtree->GetIntersecting (p, p, locels); ne = locels.Size(); } else @@ -4988,6 +5148,14 @@ namespace netgen //(*testout) << "p " << p << endl; //(*testout) << "velement " << velement << endl; + if (!GetNE() && GetNSE() ) + { + lami[0] = vlam[0]; + lami[1] = vlam[1]; + lami[2] = vlam[2]; + return velement; + } + Array faces; topology.GetElementFaces(velement,faces); @@ -6030,7 +6198,7 @@ namespace netgen void Mesh :: SetUserData(const char * id, Array & data) { if(userdata_int.Used(id)) - delete userdata_int.Get(id); + delete userdata_int[id]; Array * newdata = new Array(data); @@ -6040,10 +6208,10 @@ namespace netgen { if(userdata_int.Used(id)) { - if(data.Size() < (*userdata_int.Get(id)).Size()+shift) - data.SetSize((*userdata_int.Get(id)).Size()+shift); - for(int i=0; i<(*userdata_int.Get(id)).Size(); i++) - data[i+shift] = (*userdata_int.Get(id))[i]; + if(data.Size() < (*userdata_int[id]).Size()+shift) + data.SetSize((*userdata_int[id]).Size()+shift); + for(int i=0; i<(*userdata_int[id]).Size(); i++) + data[i+shift] = (*userdata_int[id])[i]; return true; } else @@ -6055,7 +6223,7 @@ namespace netgen void Mesh :: SetUserData(const char * id, Array & data) { if(userdata_double.Used(id)) - delete userdata_double.Get(id); + delete userdata_double[id]; Array * newdata = new Array(data); @@ -6065,10 +6233,10 @@ namespace netgen { if(userdata_double.Used(id)) { - if(data.Size() < (*userdata_double.Get(id)).Size()+shift) - data.SetSize((*userdata_double.Get(id)).Size()+shift); - for(int i=0; i<(*userdata_double.Get(id)).Size(); i++) - data[i+shift] = (*userdata_double.Get(id))[i]; + if(data.Size() < (*userdata_double[id]).Size()+shift) + data.SetSize((*userdata_double[id]).Size()+shift); + for(int i=0; i<(*userdata_double[id]).Size(); i++) + data[i+shift] = (*userdata_double[id])[i]; return true; } else