1. SIGMA
  2. Components
  3. MOAB

Commits

Iulian Grindeanu  committed 915e3f1 Merge

Merge branch 'iulian07/fix_spherical' into develop

Comments (0)

Files changed (3)

File tools/mbcoupler/ElemUtil.cpp Modified

View file
  • Ignore whitespace
  • Hide word diff
   {
     // project the vertices to the plane tangent at first vertex
     v1=vertex[0]; // member data
-    double v1v1= v1%v1;
+    double v1v1= v1%v1; // this is 1, in general, for unit sphere meshes
     for (int j=1; j<4; j++)
     {
-      CartVect vnew =v1v1/(vertex[j]%v1)*vertex[j]; // so that (vnew-v1)%v1 is 0
+      // first, bring all vertices in the gnomonic plane
+      //  the new vertex will intersect the plane at vnew
+      //   so that (vnew-v1)%v1 is 0 ( vnew is in the tangent plane, i.e. normal to v1 )
+      // pos is the old position of the vertex, and it is in general on the sphere
+      // vnew = alfa*pos;  (alfa*pos-v1)%v1 = 0  <=> alfa*(pos%v1)=v1%v1 <=> alfa = v1v1/(pos%v1)
+      //  <=> vnew = ( v1v1/(pos%v1) )*pos
+      CartVect vnew =v1v1/(vertex[j]%v1)*vertex[j];
       vertex[j]=vnew;
     }
     // will compute a transf matrix, such that a new point will be transformed with
-    // newpos =  transf * (pos-v1);
+    // newpos =  transf * (vnew-v1), and it will be a point in the 2d plane
+    // the transformation matrix will be oriented in such a way that orientation will be positive
     CartVect vx = vertex[1]-v1; // this will become Ox axis
+    // z axis will be along v1, in such a way that orientation of the quad is positive
+    // look at the first 2 edges
+    CartVect vz = vx*(vertex[2]-vertex[1]);
+    vz = vz/vz.length();
+
     vx = vx/vx.length();
-    CartVect vz = -v1/v1.length();// z will point down; with this, if the second vertex is along x,
-                                   // the det will be positive for typical HOMME meshes
-                                   // this allows to simplify the ievaluate,
-                                   // so we can just reuse Map::ievaluate
+
     CartVect vy = vz*vx;
     transf = Matrix3(vx[0], vx[1], vx[2], vy[0], vy[1], vy[2], vz[0], vz[1], vz[2]);
     vertex[0]= CartVect(0.);
       vertex[j] = transf*(vertex[j]-v1);
   }
 
-   CartVect SphericalQuad::ievaluate(double tol, const CartVect& x) const
+   CartVect SphericalQuad::ievaluate( const CartVect& x, double tol, const CartVect& x0) const
    {
-     // project to the plane tangent at first vertex
-     //CartVect v1=vertex[0];
+     // project to the plane tangent at first vertex (gnomonic projection)
      double v1v1= v1%v1;
-     CartVect vnew =v1v1/(x%v1)*x; // so that (x-v1)%v1 is 0
+     CartVect vnew =v1v1/(x%v1)*x; // so that (vnew-v1)%v1 is 0
      vnew =  transf*(vnew-v1);
      // det will be positive now
-     return Map::ievaluate(vnew, tol);
+     return Map::ievaluate(vnew, tol, x0);
    }
 
    bool SphericalQuad::inside_box(const CartVect & pos, double & tol) const

File tools/mbcoupler/ElemUtil.hpp Modified

View file
  • Ignore whitespace
  • Hide word diff
       virtual ~LinearTet();
       /* Override the evaluation routines to take advantage of the properties of P1. */
       virtual CartVect evaluate(const CartVect& xi) const {return this->vertex[0] + this->T*xi;};
-      using Map::ievaluate;
       virtual CartVect ievaluate(const CartVect& x) const {return this->T_inverse*(x-this->vertex[0]);};
       virtual Matrix3  jacobian(const CartVect& )  const {return this->T;};
       virtual Matrix3  ijacobian(const CartVect& ) const {return this->T_inverse;};
       virtual ~SpectralHex();
       void set_gl_points( double * x, double * y, double *z) ;
       virtual CartVect evaluate( const CartVect& xi ) const;
-      using Map::ievaluate;
       virtual CartVect ievaluate( double abs_eps, const CartVect& x) const;
       virtual Matrix3  jacobian(const CartVect& xi) const;
       double   evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const;
 
     };// class LinearQuad
 
-    /**\brief Shape function space for bilinear quadrilateral on sphere, obtained from the canonical linear (affine) functions. */
+    /**\brief Shape function space for bilinear quadrilateral on sphere, obtained from the
+     *  canonical linear (affine) functions.
+     *  It is mapped using gnomonic projection to a plane tangent at the first vertex
+     *  It works well for edges that are great circle arcs; RLL meshes  do not have this property, but
+     *  HOMME or MPAS meshes do have it */
     class SphericalQuad : public LinearQuad {
     public:
       SphericalQuad(const std::vector<CartVect>& vertices);
       virtual ~SphericalQuad() {};
       virtual bool inside_box(const CartVect & pos, double & tol) const;
-      using Map::ievaluate;
-      virtual CartVect ievaluate( double abs_eps, const CartVect& x) const;
+      CartVect ievaluate( const CartVect& x, double tol, const CartVect& x0 = CartVect(0.0)) const;
     protected:
       CartVect v1;
       Matrix3 transf; // so will have a lot of stuff, including the transf to a coordinate system
         virtual ~SpectralQuad();
         void set_gl_points( double * x, double * y, double *z) ;
         virtual CartVect evaluate( const CartVect& xi ) const;// a 2d, so 3rd component is 0, always
-        using Map::ievaluate;
         virtual CartVect ievaluate(const CartVect& x) const; //a 2d, so 3rd component is 0, always
         virtual Matrix3  jacobian(const CartVect& xi) const;
         double   evaluate_scalar_field(const CartVect& xi, const double *field_vertex_values) const;

File tools/mbcoupler/ElementTest.cpp Modified

View file
  • Ignore whitespace
  • Hide word diff
 void test_hex();
 void test_spectral_hex();
 void test_spectral_quad();
+void test_spherical_quad();
 
 int main()
 {
   rval += RUN_TEST(test_hex);
   rval += RUN_TEST(test_spectral_hex);
   rval += RUN_TEST(test_spectral_quad);
+  rval += RUN_TEST(test_spherical_quad);
   return rval;
 }
 
 
   delete mb;
 }
+void test_spherical_quad()
+{
+  // example from one coupler test, run like this
+  // ./mbcoupler_test -meshes sphere_16p.h5m mpas_p8.h5m -itag vertex_field -meth 4 -outfile dd.h5m
+  // method 4 is spherical
+  double positions[] =
+  {
+   -0.88882388032987436, -0.069951956448441419, 0.45287838714646161,
+   -0.88226455385461389, -0.13973697758043971, 0.4495362433757738,
+   -0.84497006020160348, -0.13383011007602069, 0.51779831884618843,
+   -0.85072691325794214, -0.066953660115039074, 0.52132612293631853
+  };
+  CartVect x(-0.85408569769999998, -0.12391301439999999, 0.50515659540000002);
+  std::vector<CartVect> vertices;
+  for (int i=0; i<4; i++)
+    vertices.push_back(CartVect(positions+3*i));
+
+  moab::Element::SphericalQuad squad(vertices);
+  double tol(0.0001);
+  if (squad.inside_box(x, tol))
+  {
+   CartVect nat_par = squad.ievaluate(x, 0.000001);
+   std::cout<< nat_par << "\n";
+  }
+  std::cout << "success...\n";
+}