Skip to content
Snippets Groups Projects
Select Git revision
  • 52ecae224ce07c05015238dbf14c5f482861526b
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

STensor3.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    GFace.cpp 15.42 KiB
    // $Id: GFace.cpp,v 1.41 2008-01-21 22:16:04 geuzaine Exp $
    //
    // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
    //
    // This program is free software; you can redistribute it and/or modify
    // it under the terms of the GNU General Public License as published by
    // the Free Software Foundation; either version 2 of the License, or
    // (at your option) any later version.
    //
    // This program is distributed in the hope that it will be useful,
    // but WITHOUT ANY WARRANTY; without even the implied warranty of
    // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    // GNU General Public License for more details.
    //
    // You should have received a copy of the GNU General Public License
    // along with this program; if not, write to the Free Software
    // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
    // USA.
    // 
    // Please report all bugs and problems to <gmsh@geuz.org>.
    
    #include "GModel.h"
    #include "GFace.h"
    #include "GEdge.h"
    #include "MElement.h"
    #include "Message.h"
    #include "Numeric.h"
    #include "Context.h"
    
    #if defined(HAVE_GSL)
    #include <gsl/gsl_vector.h>
    #include <gsl/gsl_linalg.h>
    #else
    #define NRANSI
    #include "nrutil.h"
    void dsvdcmp(double **a, int m, int n, double w[], double **v);
    #endif
    
    extern Context_T CTX;
    
    GFace::GFace(GModel *model, int tag) : GEntity(model, tag), r1(0), r2(0) 
    {
      meshStatistics.status = GFace::PENDING;
      resetMeshAttributes();
    }
    
    GFace::~GFace()
    { 
      std::list<GEdge*>::iterator it = l_edges.begin();
    
      while (it != l_edges.end()){
        (*it)->delFace(this);
        ++it;
      }
    
      for(unsigned int i = 0; i < mesh_vertices.size(); i++) 
        delete mesh_vertices[i];
    
      for(unsigned int i = 0; i < triangles.size(); i++) 
        delete triangles[i];
    
      for(unsigned int i = 0; i < quadrangles.size(); i++) 
        delete quadrangles[i];
    }
    
    void GFace::resetMeshAttributes()
    {
      meshAttributes.recombine = 0;
      meshAttributes.recombineAngle = 0.;
      meshAttributes.Method = LIBRE;
      meshAttributes.transfiniteArrangement = 0;
      meshAttributes.extrude = 0;
    }
    
    SBoundingBox3d GFace::bounds() const
    {
      SBoundingBox3d res;
      if(geomType() != DiscreteSurface){
        std::list<GEdge*>::const_iterator it = l_edges.begin();
        for(; it != l_edges.end(); it++)
          res += (*it)->bounds();
      }
      else{
        for(unsigned int i = 0; i < mesh_vertices.size(); i++)
          res += mesh_vertices[i]->point();
      }
      return res;
    }
    
    std::list<GVertex*> GFace::vertices() const
    {
      std::list<GEdge*>::const_iterator it = l_edges.begin();
      std::list<GVertex*>ret;
      while (it != l_edges.end()){
        GVertex *v1 = (*it)->getBeginVertex();
        GVertex *v2 = (*it)->getEndVertex();
        if(v1 && std::find(ret.begin(), ret.end(), v1) == ret.end())
          ret.push_back(v1);
        if(v2 && std::find(ret.begin(), ret.end(), v2) == ret.end())
          ret.push_back(v2);
        ++it;
      }
      return ret;
    }
    
    void GFace::setVisibility(char val, bool recursive)
    {
      GEntity::setVisibility(val);
      if(recursive){
        std::list<GEdge*>::iterator it = l_edges.begin();
        while (it != l_edges.end()){
          (*it)->setVisibility(val, recursive);
          ++it;
        }
      }
    }
    
    void GFace::recomputeMeshPartitions()
    {
      for(unsigned int i = 0; i < triangles.size(); i++) {
        int part = triangles[i]->getPartition();
        if(part) model()->getMeshPartitions().insert(part);
      }
      for(unsigned int i = 0; i < quadrangles.size(); i++) {
        int part = quadrangles[i]->getPartition();
        if(part) model()->getMeshPartitions().insert(part);
      }
    }
    
    void GFace::deleteMeshPartitions()
    {
      for(unsigned int i = 0; i < triangles.size(); i++)
        triangles[i]->setPartition(0);
      for(unsigned int i = 0; i < quadrangles.size(); i++)
        quadrangles[i]->setPartition(0);
    }
    
    std::string GFace::getAdditionalInfoString()
    {
      if(l_edges.empty()) return std::string("");
    
      char tmp[256];
      if(l_edges.size() > 10){
        sprintf(tmp, "{%d, ..., %d}", (*l_edges.begin())->tag(), (*l_edges.end())->tag());
        return std::string(tmp);
      }
    
      std::string str("");
      std::list<GEdge*>::const_iterator it = l_edges.begin();
      str += "{";
      for(; it != l_edges.end(); it++){
        if(it != l_edges.begin()) str += ",";
        sprintf(tmp, "%d", (*it)->tag());
        str += tmp;
      }
      str += "}";
      return str;
    }
    
    void GFace::computeMeanPlane()
    {
      std::vector<SPoint3> pts;
      std::list<GVertex*> verts = vertices();
      std::list<GVertex*>::const_iterator itv = verts.begin();
      for(; itv != verts.end(); itv++){
        const GVertex *v = *itv; 
        pts.push_back(SPoint3(v->x(), v->y(), v->z()));
      }
    
      if(pts.size() < 3){
        Msg(INFO, "Adding edge points to compute mean plane of face %d", tag());
        std::list<GEdge*> edg = edges();
        std::list<GEdge*>::const_iterator ite = edg.begin();
        for(; ite != edg.end(); ite++){
          const GEdge *e = *ite; 
          Range<double> b = e->parBounds(0);
          GPoint p1 = e->point(b.low() + 0.333 * (b.high() - b.low()));
          pts.push_back(SPoint3(p1.x(), p1.y(), p1.z()));
          GPoint p2 = e->point(b.low() + 0.666 * (b.high() - b.low()));
          pts.push_back(SPoint3(p2.x(), p2.y(), p2.z()));
        }
      }
    
      computeMeanPlane(pts);
    }
    
    void GFace::computeMeanPlane(const std::vector<MVertex*> &points)
    {
      std::vector<SPoint3> pts;
      for(unsigned int i = 0; i < points.size(); i++)
        pts.push_back(SPoint3(points[i]->x(), points[i]->y(), points[i]->z()));
      computeMeanPlane(pts);
    }
    
    void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
    {
      // The concept of a mean plane computed in the sense of least
      // squares is fine for plane surfaces(!), but not really the best
      // one for non-plane surfaces. Indeed, imagine a quarter of a circle
      // extruded along a very small height: the mean plane will be in the
      // plane of the circle... Until we implement something better, there
      // is a test in this routine that checks the coherence of the
      // computation for non-plane surfaces and reverts to a basic mean
      // plane approximatation if the check fails.
    
      double xm = 0., ym = 0., zm = 0.;
      int ndata = points.size();
      int na = 3;
      for(int i = 0; i < ndata; i++) {
        xm += points[i].x();
        ym += points[i].y();
        zm += points[i].z();
      }
      xm /= (double)ndata;
      ym /= (double)ndata;
      zm /= (double)ndata;
    
      int min;
      double res[4], svd[3];
    #if defined(HAVE_GSL)
      gsl_matrix *U = gsl_matrix_alloc(ndata, na);
      gsl_matrix *V = gsl_matrix_alloc(na, na);
      gsl_vector *W = gsl_vector_alloc(na);
      gsl_vector *TMPVEC = gsl_vector_alloc(na);
      for(int i = 0; i < ndata; i++) {
        gsl_matrix_set(U, i, 0, points[i].x() - xm);
        gsl_matrix_set(U, i, 1, points[i].y() - ym);
        gsl_matrix_set(U, i, 2, points[i].z() - zm);
      }
      gsl_linalg_SV_decomp(U, V, W, TMPVEC);
      svd[0] = gsl_vector_get(W, 0);
      svd[1] = gsl_vector_get(W, 1);
      svd[2] = gsl_vector_get(W, 2);
      if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
        min = 0;
      else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
        min = 1;
      else
        min = 2;
      res[0] = gsl_matrix_get(V, 0, min);
      res[1] = gsl_matrix_get(V, 1, min);
      res[2] = gsl_matrix_get(V, 2, min);
      norme(res);
      gsl_matrix_free(U);
      gsl_matrix_free(V);
      gsl_vector_free(W);
      gsl_vector_free(TMPVEC);
    #else
      double **U = dmatrix(1, ndata, 1, na);
      double **V = dmatrix(1, na, 1, na);
      double *W = dvector(1, na);
      for(int i = 0; i < ndata; i++) {
        U[i + 1][1] = points[i].x() - xm;
        U[i + 1][2] = points[i].y() - ym;
        U[i + 1][3] = points[i].z() - zm;
      }
      dsvdcmp(U, ndata, na, W, V);
      if(fabs(W[1]) < fabs(W[2]) && fabs(W[1]) < fabs(W[3]))
        min = 1;
      else if(fabs(W[2]) < fabs(W[1]) && fabs(W[2]) < fabs(W[3]))
        min = 2;
      else
        min = 3;
      svd[0] = W[1];
      svd[1] = W[2];
      svd[2] = W[3];
      res[0] = V[1][min];
      res[1] = V[2][min];
      res[2] = V[3][min];
      norme(res);
      free_dmatrix(U, 1, ndata, 1, na);
      free_dmatrix(V, 1, na, 1, na);
      free_dvector(W, 1, na);
    #endif
    
      double ex[3], t1[3], t2[3];
    
      // check coherence of results for non-plane surfaces
      if(geomType() != Plane && geomType() != DiscreteSurface) {
        double res2[3], c[3], cosc, sinc, angplan;
        double eps = 1.e-3;
    
        GPoint v1 = point(0.5, 0.5);
        GPoint v2 = point(0.5 + eps, 0.5);
        GPoint v3 = point(0.5, 0.5 + eps);
        t1[0] = v2.x() - v1.x();
        t1[1] = v2.y() - v1.y();
        t1[2] = v2.z() - v1.z();
        t2[0] = v3.x() - v1.x();
        t2[1] = v3.y() - v1.y();
        t2[2] = v3.z() - v1.z();
        norme(t1);
        norme(t2);
        // prodve(t1, t2, res2);
        // Warning: the rest of the code assumes res = t2 x t1, not t1 x t2 (WTF?)
        prodve(t2, t1, res2); 
        norme(res2);
        prodve(res, res2, c);
        prosca(res, res2, &cosc);
        sinc = sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
        angplan = myatan2(sinc, cosc);
        angplan = angle_02pi(angplan) * 180. / Pi;
        if((angplan > 70 && angplan < 110) || (angplan > 260 && angplan < 280)) {
          Msg(INFO, "SVD failed (angle=%g): using rough algo...", angplan);
          res[0] = res2[0];
          res[1] = res2[1];
          res[2] = res2[2];
          goto end;
        }
      }
    
      ex[0] = ex[1] = ex[2] = 0.0;
      if(res[0] == 0.)
        ex[0] = 1.0;
      else if(res[1] == 0.)
        ex[1] = 1.0;
      else
        ex[2] = 1.0;
    
      prodve(res, ex, t1);
      norme(t1);
      prodve(t1, res, t2);
      norme(t2);
    
    end:
      res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
    
      for(int i = 0; i < 3; i++)
        meanPlane.plan[0][i] = t1[i];
      for(int i = 0; i < 3; i++)
        meanPlane.plan[1][i] = t2[i];
      for(int i = 0; i < 3; i++)
        meanPlane.plan[2][i] = res[i];
    
      meanPlane.a = res[0];
      meanPlane.b = res[1];
      meanPlane.c = res[2];
      meanPlane.d = res[3];
    
      meanPlane.x = meanPlane.y = meanPlane.z = 0.;
      if(fabs(meanPlane.a) >= fabs(meanPlane.b) && 
         fabs(meanPlane.a) >= fabs(meanPlane.c) ){
        meanPlane.x = meanPlane.d / meanPlane.a;
      }
      else if(fabs(meanPlane.b) >= fabs(meanPlane.a) && 
    	  fabs(meanPlane.b) >= fabs(meanPlane.c)){
        meanPlane.y = meanPlane.d / meanPlane.b;
      }
      else{
        meanPlane.z = meanPlane.d / meanPlane.c;
      }
      
      Msg(DEBUG1, "Surface: %d", tag());
      Msg(DEBUG2, "SVD    : %g,%g,%g (min=%d)", svd[0], svd[1], svd[2], min);
      Msg(DEBUG2, "Plane  : (%g x + %g y + %g z = %g)", 
          meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d);
      Msg(DEBUG2, "Normal : (%g , %g , %g )", 
          meanPlane.a, meanPlane.b, meanPlane.c);
      Msg(DEBUG3, "t1     : (%g , %g , %g )", t1[0], t1[1], t1[2]);
      Msg(DEBUG3, "t2     : (%g , %g , %g )", t2[0], t2[1], t2[2]);
      Msg(DEBUG3, "pt     : (%g , %g , %g )", 
          meanPlane.x, meanPlane.y, meanPlane.z);
    
      //check coherence for plane surfaces
      if(geomType() == Plane) {
        SBoundingBox3d bb = bounds();
        double lc = norm(SVector3(bb.max(), bb.min()));
        std::list<GVertex*> verts = vertices();
        std::list<GVertex*>::const_iterator itv = verts.begin();
        for(; itv != verts.end(); itv++){
          const GVertex *v = *itv; 
          double d = meanPlane.a * v->x() + meanPlane.b * v->y() + 
    	meanPlane.c * v->z() - meanPlane.d;
          if(fabs(d) > lc * 1.e-3) {
    	Msg(GERROR1, "Plane surface %d (%gx+%gy+%gz+%g=0) is not plane!",
    	    tag(), meanPlane.a, meanPlane.b, meanPlane.c, meanPlane.d);
    	Msg(GERROR3, "Control point %d = (%g,%g,%g), val=%g",
    	    v->tag(), v->x(), v->y(), v->z(), d);
    	return;
          }
        }
      }
    }
    
    void GFace::getMeanPlaneData(double VX[3], double VY[3], 
    			     double &x, double &y, double &z) const
    {
      VX[0] = meanPlane.plan[0][0];
      VX[1] = meanPlane.plan[0][1];
      VX[2] = meanPlane.plan[0][2];
      VY[0] = meanPlane.plan[1][0];
      VY[1] = meanPlane.plan[1][1];
      VY[2] = meanPlane.plan[1][2];
      x = meanPlane.x;  
      y = meanPlane.y;  
      z = meanPlane.z;  
    }
    
    void GFace::getMeanPlaneData(double plan[3][3]) const
    {
      for(int i = 0; i < 3; i++)
        for(int j = 0; j < 3; j++)
          plan[i][j] = meanPlane.plan[i][j];
    }
    
    double GFace::curvature (const SPoint2 &param) const
    {
      if (geomType() == Plane) return 0;
    
      // X=X(u,v) Y=Y(u,v) Z=Z(u,v)
      // curv = div n = dnx/dx + dny/dy + dnz/dz
      // dnx/dx = dnx/du du/dx + dnx/dv dv/dx
    
      const double eps = 1.e-3;
    
      Pair<SVector3,SVector3> der = firstDer(param);
    
      SVector3 du = der.first();
      SVector3 dv = der.second();
      SVector3 nml = crossprod(du, dv);
    
      double detJ = norm(nml);
    
      du.normalize();
      dv.normalize();
    
      SVector3 n1 = normal(SPoint2(param.x() - eps, param.y()));
      SVector3 n2 = normal(SPoint2(param.x() + eps, param.y()));
      SVector3 n3 = normal(SPoint2(param.x(), param.y() - eps));
      SVector3 n4 = normal(SPoint2(param.x(), param.y() + eps));
    
      SVector3 dndu = 500 * (n2 - n1);
      SVector3 dndv = 500 * (n4 - n3);
    
      double c = fabs(dot(dndu, du) +  dot(dndv, dv)) / detJ; 
    
      // Msg(INFO, "c = %g detJ %g", c, detJ);
    
      return c;
    }
    
    void GFace::XYZtoUV(const double X, const double Y, const double Z, 
    		    double &U, double &V, const double relax,
    		    const bool onSurface) const
    {
      const double Precision = 1.e-8;
      const int MaxIter = 25;
      const int NumInitGuess = 11;
    
      double Unew = 0., Vnew = 0., err, err2;
      int iter;
      double mat[3][3], jac[3][3];
      double umin, umax, vmin, vmax;
      double initu[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 1., 0., 0., 1.};
      double initv[NumInitGuess] = {0.487, 0.6, 0.4, 0.7, 0.3, 0.8, 0.2, 0., 1., 0., 1.};
      
      Range<double> ru = parBounds(0);
      Range<double> rv = parBounds(1);
      umin = ru.low();
      umax = ru.high();
      vmin = rv.low();
      vmax = rv.high();
    
      for(int i = 0; i < NumInitGuess; i++){
        for(int j = 0; j < NumInitGuess; j++){
        
          U = initu[i];
          V = initv[j];
          err = 1.0;
          iter = 1;
    
          GPoint P = point(U,V);
          err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
          if (err2 < 1.e-8 * CTX.lc)return;
    
    
          while(err > Precision && iter < MaxIter) {
    	P = point(U,V);
    	Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V)); 
    	mat[0][0] = der.left().x();
    	mat[0][1] = der.left().y();
    	mat[0][2] = der.left().z();
    	mat[1][0] = der.right().x();
    	mat[1][1] = der.right().y();
    	mat[1][2] = der.right().z();
    	mat[2][0] = 0.;
    	mat[2][1] = 0.;
    	mat[2][2] = 0.;
    	invert_singular_matrix3x3(mat, jac);
    	
    	Unew = U + relax *
    	  (jac[0][0] * (X - P.x()) + jac[1][0] * (Y - P.y()) +
    	   jac[2][0] * (Z - P.z()));
    	Vnew = V + relax * 
    	  (jac[0][1] * (X - P.x()) + jac[1][1] * (Y - P.y()) +
    	   jac[2][1] * (Z - P.z()));
    	
    	err = DSQR(Unew - U) + DSQR(Vnew - V);
    	err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
    	iter++;
    	U = Unew;
    	V = Vnew;
          }
          
          if(iter < MaxIter && err <= Precision && 
    	 Unew <= umax && Vnew <= vmax && 
    	 Unew >= umin && Vnew >= vmin){
    	if (onSurface && err2 > 1.e-4 * CTX.lc)
    	  Msg(WARNING,"Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g", 
    	      i, j, err, iter, err2);
    	return;	
          }
        }
      }
      
      if(relax < 1.e-6)
        Msg(GERROR, "Could not converge: surface mesh will be wrong");
      else {
        Msg(INFO, "point %g %g %g : Relaxation factor = %g",X,Y,Z, 0.75 * relax);
        XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
      }  
    }
    
    SPoint2 GFace::parFromPoint(const SPoint3 &p) const
    {
      double U,V;
      XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
      return SPoint2(U,V);
    }
    
    void GFace::computeGraphicsRep(int nu, int nv)
    {
      _graphicsRep.resize(nu);
      for(int i = 0; i < nu; i++) _graphicsRep[i].resize(nv);
    
      Range<double> ubounds = parBounds(0);
      Range<double> vbounds = parBounds(1);
      double umin = ubounds.low(), umax = ubounds.high();
      double vmin = vbounds.low(), vmax = vbounds.high();
    
      for(int i = 0; i < nu; i++){
        for(int j = 0; j < nv; j++){
          double u = umin + (double)i/(double)(nu-1) * (umax - umin);
          double v = vmin + (double)j/(double)(nv-1) * (vmax - vmin);
          struct graphics_point gp;
          GPoint p = point(u, v);
          gp.xyz[0] = p.x();
          gp.xyz[1] = p.y();
          gp.xyz[2] = p.z();
          SVector3 n = normal(SPoint2(u, v));
          gp.n[0] = n.x();
          gp.n[1] = n.y();
          gp.n[2] = n.z();
          _graphicsRep[i][j] = gp;
        }
      }
    }