# Commits

committed 915e3f1 Merge

Merge branch 'iulian07/fix_spherical' into develop

# File tools/mbcoupler/ElemUtil.cpp Modified

• 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

• 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

• 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";`
`+}`