[MOAB-dev] r2836 - in MOAB/trunk: . tools/mbcoupler
kraftche at cae.wisc.edu
kraftche at cae.wisc.edu
Fri Apr 17 14:44:03 CDT 2009
Author: kraftche
Date: 2009-04-17 14:44:02 -0500 (Fri, 17 Apr 2009)
New Revision: 2836
Added:
MOAB/trunk/MBElemUtil.cpp
MOAB/trunk/MBElemUtil.hpp
Removed:
MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
Modified:
MOAB/trunk/MBGeomUtil.hpp
MOAB/trunk/Makefile.am
MOAB/trunk/tools/mbcoupler/Makefile.am
Log:
move MBElemUtil back into MOAB proper
Copied: MOAB/trunk/MBElemUtil.cpp (from rev 2835, MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp)
===================================================================
--- MOAB/trunk/MBElemUtil.cpp (rev 0)
+++ MOAB/trunk/MBElemUtil.cpp 2009-04-17 19:44:02 UTC (rev 2836)
@@ -0,0 +1,297 @@
+#include <iostream>
+#include <limits>
+#include <assert.h>
+
+#include "MBElemUtil.hpp"
+#include "MBMatrix3.hpp"
+#include "types.h"
+
+namespace MBElemUtil {
+
+/**\brief Class representing a 3-D mapping function (e.g. shape function for volume element) */
+class VolMap {
+ public:
+ /**\brief Return $\vec \xi$ corresponding to logical center of element */
+ virtual MBCartVect center_xi() const = 0;
+ /**\brief Evaluate mapping function (calculate $\vec x = F($\vec \xi)$ )*/
+ virtual MBCartVect evaluate( const MBCartVect& xi ) const = 0;
+ /**\brief Evaluate Jacobian of mapping function */
+ virtual MBMatrix3 jacobian( const MBCartVect& xi ) const = 0;
+ /**\brief Evaluate inverse of mapping function (calculate $\vec \xi = F^-1($\vec x)$ )*/
+ bool solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const ;
+};
+
+bool VolMap::solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const
+{
+ const double error_tol_sqr = tol*tol;
+ double det;
+ xi = center_xi();
+ MBCartVect delta = evaluate(xi) - x;
+ MBMatrix3 J;
+ while (delta % delta > error_tol_sqr) {
+ J = jacobian(xi);
+ det = J.determinant();
+ if (det < std::numeric_limits<double>::epsilon())
+ return false;
+ xi -= J.inverse(1.0/det) * delta;
+ delta = evaluate( xi ) - x;
+ }
+ return true;
+}
+
+/**\brief Shape function for trilinear hexahedron */
+class LinearHexMap : public VolMap {
+ public:
+ LinearHexMap( const MBCartVect* corner_coords ) : corners(corner_coords) {}
+ virtual MBCartVect center_xi() const;
+ virtual MBCartVect evaluate( const MBCartVect& xi ) const;
+ virtual MBMatrix3 jacobian( const MBCartVect& xi ) const;
+ private:
+ const MBCartVect* corners;
+ static const double corner_xi[8][3];
+};
+
+const double LinearHexMap::corner_xi[8][3] = { { -1, -1, -1 },
+ { 1, -1, -1 },
+ { 1, 1, -1 },
+ { -1, 1, -1 },
+ { -1, -1, 1 },
+ { 1, -1, 1 },
+ { 1, 1, 1 },
+ { -1, 1, 1 } };
+MBCartVect LinearHexMap::center_xi() const
+ { return MBCartVect(0.0); }
+
+MBCartVect LinearHexMap::evaluate( const MBCartVect& xi ) const
+{
+ MBCartVect x(0.0);
+ for (unsigned i = 0; i < 8; ++i) {
+ const double N_i = (1 + xi[0]*corner_xi[i][0])
+ * (1 + xi[1]*corner_xi[i][1])
+ * (1 + xi[2]*corner_xi[i][2]);
+ x += N_i * corners[i];
+ }
+ x *= 0.125;
+ return x;
+}
+
+MBMatrix3 LinearHexMap::jacobian( const MBCartVect& xi ) const
+{
+ MBMatrix3 J(0.0);
+ for (unsigned i = 0; i < 8; ++i) {
+ const double xi_p = 1 + xi[0]*corner_xi[i][0];
+ const double eta_p = 1 + xi[1]*corner_xi[i][1];
+ const double zeta_p = 1 + xi[2]*corner_xi[i][2];
+ const double dNi_dxi = corner_xi[i][0] * eta_p * zeta_p;
+ const double dNi_deta = corner_xi[i][1] * xi_p * zeta_p;
+ const double dNi_dzeta = corner_xi[i][2] * xi_p * eta_p;
+ J(0,0) += dNi_dxi * corners[i][0];
+ J(1,0) += dNi_dxi * corners[i][1];
+ J(2,0) += dNi_dxi * corners[i][2];
+ J(0,1) += dNi_deta * corners[i][0];
+ J(1,1) += dNi_deta * corners[i][1];
+ J(2,1) += dNi_deta * corners[i][2];
+ J(0,2) += dNi_dzeta * corners[i][0];
+ J(1,2) += dNi_dzeta * corners[i][1];
+ J(2,2) += dNi_dzeta * corners[i][2];
+ }
+ return J *= 0.125;
+}
+
+bool nat_coords_trilinear_hex( const MBCartVect* corner_coords,
+ const MBCartVect& x,
+ MBCartVect& xi,
+ double tol )
+{
+ return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
+}
+
+
+//
+// nat_coords_trilinear_hex2
+// Duplicate functionality of nat_coords_trilinear_hex using hex_findpt
+//
+void nat_coords_trilinear_hex2(const MBCartVect hex[8],
+ const MBCartVect& xyz,
+ MBCartVect &ncoords,
+ double etol)
+
+{
+ const int ndim = 3;
+ const int nverts = 8;
+ const int vertMap[nverts] = {0,1,3,2, 4,5,7,6}; //Map from nat to lex ordering
+
+ const int n = 2; //linear
+ real coords[ndim*nverts]; //buffer
+
+ real *xm[ndim];
+ for(int i=0; i<ndim; i++)
+ xm[i] = coords + i*nverts;
+
+ //stuff hex into coords
+ for(int i=0; i<nverts; i++){
+ real vcoord[ndim];
+ hex[i].get(vcoord);
+
+ for(int d=0; d<ndim; d++)
+ coords[d*nverts + vertMap[i]] = vcoord[d];
+
+ }
+
+ double dist = 0.0;
+ MBElemUtil::hex_findpt(xm, n, xyz, ncoords, dist);
+ if (3*EPS < dist) {
+ // outside element, set extremal values to something outside range
+ for (int j = 0; j < 3; j++) {
+ if (ncoords[j] < (-1.0-etol) || ncoords[j] > (1.0+etol))
+ ncoords[j] *= 10;
+ }
+ }
+
+}
+
+bool point_in_trilinear_hex(const MBCartVect *hex,
+ const MBCartVect& xyz,
+ double etol)
+{
+ MBCartVect xi;
+ return nat_coords_trilinear_hex( hex, xyz, xi, etol )
+ && fabs(xi[0])-1 < etol
+ && fabs(xi[1])-1 < etol
+ && fabs(xi[2])-1 < etol;
+}
+
+
+bool point_in_trilinear_hex(const MBCartVect *hex,
+ const MBCartVect& xyz,
+ const MBCartVect& box_min,
+ const MBCartVect& box_max,
+ double etol)
+{
+ // all values scaled by 2 (eliminates 3 flops)
+ const MBCartVect mid = box_max + box_min;
+ const MBCartVect dim = box_max - box_min;
+ const MBCartVect pt = 2*xyz - mid;
+ return fabs(pt[0]) - dim[0] < etol &&
+ fabs(pt[1]) - dim[1] < etol &&
+ fabs(pt[2]) - dim[2] < etol &&
+ point_in_trilinear_hex( hex, xyz, etol );
+}
+
+
+
+// Wrapper to James Lottes' findpt routines
+// hex_findpt
+// Find the parametric coordinates of an xyz point inside
+// a 3d hex spectral element with n nodes per dimension
+// xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
+// n: number of nodes per dimension -- n=2 for a linear element
+// xyz: input, point to find
+// rst: output: parametric coords of xyz inside the element. If xyz is outside the element, rst will be the coords of the closest point
+// dist: output: distance between xyz and the point with parametric coords rst
+extern "C"{
+#include "types.h"
+#include "poly.h"
+#include "tensor.h"
+#include "findpt.h"
+#include "extrafindpt.h"
+#include "errmem.h"
+}
+
+void hex_findpt(real *xm[3],
+ int n,
+ MBCartVect xyz,
+ MBCartVect &rst,
+ double &dist)
+{
+
+ //compute stuff that only depends on the order -- could be cached
+ real *z[3];
+ lagrange_data ld[3];
+ opt_data_3 data;
+
+ //triplicates
+ for(int d=0; d<3; d++){
+ z[d] = tmalloc(real, n);
+ lobatto_nodes(z[d], n);
+ lagrange_setup(&ld[d], z[d], n);
+ }
+
+ opt_alloc_3(&data, ld);
+
+ //find nearest point
+ real x_star[3];
+ xyz.get(x_star);
+
+ real r[3] = {0, 0, 0 }; // initial guess for parametric coords
+ unsigned c = opt_no_constraints_3;
+ dist = opt_findpt_3(&data, (const real **)xm, x_star, r, &c);
+ //c tells us if we landed inside the element or exactly on a face, edge, or node
+
+ //copy parametric coords back
+ rst = r;
+
+ //Clean-up (move to destructor if we decide to cache)
+ opt_free_3(&data);
+ for(int d=0; d<3; ++d)
+ lagrange_free(&ld[d]);
+ for(int d=0; d<3; ++d)
+ free(z[d]);
+}
+
+
+
+
+// hex_eval
+// Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given parametric coordinates
+// field: field values for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
+// n: number of nodes per dimension -- n=2 for a linear element
+// rst: input: parametric coords of the point where we want to evaluate the field
+// value: output: value of field at rst
+
+void hex_eval(real *field,
+ int n,
+ MBCartVect rstCartVec,
+ double &value)
+{
+ int d;
+ real rst[3];
+ rstCartVec.get(rst);
+
+ //can cache stuff below
+ lagrange_data ld[3];
+ real *z[3];
+ for(d=0;d<3;++d){
+ z[d] = tmalloc(real, n);
+ lobatto_nodes(z[d], n);
+ lagrange_setup(&ld[d], z[d], n);
+ }
+
+ //cut and paste -- see findpt.c
+ const unsigned
+ nf = n*n,
+ ne = n,
+ nw = 2*n*n + 3*n;
+ real *od_work = tmalloc(real, 6*nf + 9*ne + nw);
+
+ //piece that we shouldn't want to cache
+ for(d=0; d<3; d++){
+ lagrange_0(&ld[d], rst[d]);
+ }
+
+ value = tensor_i3(ld[0].J,ld[0].n,
+ ld[1].J,ld[1].n,
+ ld[2].J,ld[2].n,
+ field,
+ od_work);
+
+ //all this could be cached
+ for(d=0; d<3; d++){
+ free(z[d]);
+ lagrange_free(&ld[d]);
+ }
+ free(od_work);
+}
+
+} // namespace MBElemUtil
+
Property changes on: MOAB/trunk/MBElemUtil.cpp
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: MOAB/trunk/MBElemUtil.hpp (from rev 2835, MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp)
===================================================================
--- MOAB/trunk/MBElemUtil.hpp (rev 0)
+++ MOAB/trunk/MBElemUtil.hpp 2009-04-17 19:44:02 UTC (rev 2836)
@@ -0,0 +1,38 @@
+#include "MBCore.hpp"
+#include "MBCartVect.hpp"
+
+namespace MBElemUtil {
+
+ bool nat_coords_trilinear_hex(const MBCartVect* hex_corners,
+ const MBCartVect& x,
+ MBCartVect& xi,
+ double tol);
+ bool point_in_trilinear_hex(const MBCartVect *hex_corners,
+ const MBCartVect& xyz,
+ double etol);
+
+ bool point_in_trilinear_hex(const MBCartVect *hex_corners,
+ const MBCartVect& xyz,
+ const MBCartVect& box_min,
+ const MBCartVect& box_max,
+ double etol);
+
+ //wrapper to hex_findpt
+ void nat_coords_trilinear_hex2(const MBCartVect* hex_corners,
+ const MBCartVect& x,
+ MBCartVect& xi,
+ double til);
+
+
+
+ void hex_findpt(double *xm[3],
+ int n,
+ MBCartVect xyz,
+ MBCartVect& rst,
+ double& dist);
+
+ void hex_eval(double *field,
+ int n,
+ MBCartVect rst,
+ double &value);
+} // namespace MBElemUtil
Property changes on: MOAB/trunk/MBElemUtil.hpp
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: MOAB/trunk/MBGeomUtil.hpp
===================================================================
--- MOAB/trunk/MBGeomUtil.hpp 2009-04-17 16:49:56 UTC (rev 2835)
+++ MOAB/trunk/MBGeomUtil.hpp 2009-04-17 19:44:02 UTC (rev 2836)
@@ -254,30 +254,6 @@
const MBCartVect& box_center,
const MBCartVect& box_dims);
-//
-// point_in_trilinear_hex
-// Tests if a point in xyz space is within a hex element defined with
-// its eight vertex points forming a trilinear basis function. Computes
-// the natural coordinates with respect to the hex of the xyz point
-// and checks if each are between +/-1. If anyone is outside the range
-// the function returns false, otherwise it returns true.
-//
-bool point_in_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- double etol);
-
-//
-// point_in_trilinear_hex
-// Tests if a point in xyz space is within a hex element defined with
-// its eight vertex points forming a trilinear basis function. Like the
-// test above except that it gets a bounding box as arguments to filter
-// the test for acceleration.
-//
-bool point_in_trilinear_hex(MBCartVect *hex,
- MBCartVect xyz,
- MBCartVect box_min, MBCartVect box_max,
- double etol);
-
} // namespace MBGeoemtry
#endif
Modified: MOAB/trunk/Makefile.am
===================================================================
--- MOAB/trunk/Makefile.am 2009-04-17 16:49:56 UTC (rev 2835)
+++ MOAB/trunk/Makefile.am 2009-04-17 19:44:02 UTC (rev 2836)
@@ -139,6 +139,7 @@
MBCN.cpp \
MBCNArrays.hpp \
MBCartVect.cpp \
+ MBElemUtil.cpp \
MBMatrix3.cpp \
MBMatrix3.hpp \
MBCore.cpp \
@@ -245,6 +246,7 @@
MBCN_protos.h \
MBCartVect.hpp \
MBCore.hpp \
+ MBElemUtil.hpp \
MBEntityType.h \
MBEntityHandle.h \
MBError.hpp \
Deleted: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2009-04-17 16:49:56 UTC (rev 2835)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2009-04-17 19:44:02 UTC (rev 2836)
@@ -1,297 +0,0 @@
-#include <iostream>
-#include <limits>
-#include <assert.h>
-
-#include "MBElemUtil.hpp"
-#include "MBMatrix3.hpp"
-#include "types.h"
-
-namespace MBElemUtil {
-
-/**\brief Class representing a 3-D mapping function (e.g. shape function for volume element) */
-class VolMap {
- public:
- /**\brief Return $\vec \xi$ corresponding to logical center of element */
- virtual MBCartVect center_xi() const = 0;
- /**\brief Evaluate mapping function (calculate $\vec x = F($\vec \xi)$ )*/
- virtual MBCartVect evaluate( const MBCartVect& xi ) const = 0;
- /**\brief Evaluate Jacobian of mapping function */
- virtual MBMatrix3 jacobian( const MBCartVect& xi ) const = 0;
- /**\brief Evaluate inverse of mapping function (calculate $\vec \xi = F^-1($\vec x)$ )*/
- bool solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const ;
-};
-
-bool VolMap::solve_inverse( const MBCartVect& x, MBCartVect& xi, double tol ) const
-{
- const double error_tol_sqr = tol*tol;
- double det;
- xi = center_xi();
- MBCartVect delta = evaluate(xi) - x;
- MBMatrix3 J;
- while (delta % delta > error_tol_sqr) {
- J = jacobian(xi);
- det = J.determinant();
- if (det < std::numeric_limits<double>::epsilon())
- return false;
- xi -= J.inverse(1.0/det) * delta;
- delta = evaluate( xi ) - x;
- }
- return true;
-}
-
-/**\brief Shape function for trilinear hexahedron */
-class LinearHexMap : public VolMap {
- public:
- LinearHexMap( const MBCartVect* corner_coords ) : corners(corner_coords) {}
- virtual MBCartVect center_xi() const;
- virtual MBCartVect evaluate( const MBCartVect& xi ) const;
- virtual MBMatrix3 jacobian( const MBCartVect& xi ) const;
- private:
- const MBCartVect* corners;
- static const double corner_xi[8][3];
-};
-
-const double LinearHexMap::corner_xi[8][3] = { { -1, -1, -1 },
- { 1, -1, -1 },
- { 1, 1, -1 },
- { -1, 1, -1 },
- { -1, -1, 1 },
- { 1, -1, 1 },
- { 1, 1, 1 },
- { -1, 1, 1 } };
-MBCartVect LinearHexMap::center_xi() const
- { return MBCartVect(0.0); }
-
-MBCartVect LinearHexMap::evaluate( const MBCartVect& xi ) const
-{
- MBCartVect x(0.0);
- for (unsigned i = 0; i < 8; ++i) {
- const double N_i = (1 + xi[0]*corner_xi[i][0])
- * (1 + xi[1]*corner_xi[i][1])
- * (1 + xi[2]*corner_xi[i][2]);
- x += N_i * corners[i];
- }
- x *= 0.125;
- return x;
-}
-
-MBMatrix3 LinearHexMap::jacobian( const MBCartVect& xi ) const
-{
- MBMatrix3 J(0.0);
- for (unsigned i = 0; i < 8; ++i) {
- const double xi_p = 1 + xi[0]*corner_xi[i][0];
- const double eta_p = 1 + xi[1]*corner_xi[i][1];
- const double zeta_p = 1 + xi[2]*corner_xi[i][2];
- const double dNi_dxi = corner_xi[i][0] * eta_p * zeta_p;
- const double dNi_deta = corner_xi[i][1] * xi_p * zeta_p;
- const double dNi_dzeta = corner_xi[i][2] * xi_p * eta_p;
- J(0,0) += dNi_dxi * corners[i][0];
- J(1,0) += dNi_dxi * corners[i][1];
- J(2,0) += dNi_dxi * corners[i][2];
- J(0,1) += dNi_deta * corners[i][0];
- J(1,1) += dNi_deta * corners[i][1];
- J(2,1) += dNi_deta * corners[i][2];
- J(0,2) += dNi_dzeta * corners[i][0];
- J(1,2) += dNi_dzeta * corners[i][1];
- J(2,2) += dNi_dzeta * corners[i][2];
- }
- return J *= 0.125;
-}
-
-bool nat_coords_trilinear_hex( const MBCartVect* corner_coords,
- const MBCartVect& x,
- MBCartVect& xi,
- double tol )
-{
- return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
-}
-
-
-//
-// nat_coords_trilinear_hex2
-// Duplicate functionality of nat_coords_trilinear_hex using hex_findpt
-//
-void nat_coords_trilinear_hex2(const MBCartVect hex[8],
- const MBCartVect& xyz,
- MBCartVect &ncoords,
- double etol)
-
-{
- const int ndim = 3;
- const int nverts = 8;
- const int vertMap[nverts] = {0,1,3,2, 4,5,7,6}; //Map from nat to lex ordering
-
- const int n = 2; //linear
- real coords[ndim*nverts]; //buffer
-
- real *xm[ndim];
- for(int i=0; i<ndim; i++)
- xm[i] = coords + i*nverts;
-
- //stuff hex into coords
- for(int i=0; i<nverts; i++){
- real vcoord[ndim];
- hex[i].get(vcoord);
-
- for(int d=0; d<ndim; d++)
- coords[d*nverts + vertMap[i]] = vcoord[d];
-
- }
-
- double dist = 0.0;
- MBElemUtil::hex_findpt(xm, n, xyz, ncoords, dist);
- if (3*EPS < dist) {
- // outside element, set extremal values to something outside range
- for (int j = 0; j < 3; j++) {
- if (ncoords[j] < (-1.0-etol) || ncoords[j] > (1.0+etol))
- ncoords[j] *= 10;
- }
- }
-
-}
-
-bool point_in_trilinear_hex(const MBCartVect *hex,
- const MBCartVect& xyz,
- double etol)
-{
- MBCartVect xi;
- return nat_coords_trilinear_hex( hex, xyz, xi, etol )
- && fabs(xi[0])-1 < etol
- && fabs(xi[1])-1 < etol
- && fabs(xi[2])-1 < etol;
-}
-
-
-bool point_in_trilinear_hex(const MBCartVect *hex,
- const MBCartVect& xyz,
- const MBCartVect& box_min,
- const MBCartVect& box_max,
- double etol)
-{
- // all values scaled by 2 (eliminates 3 flops)
- const MBCartVect mid = box_max + box_min;
- const MBCartVect dim = box_max - box_min;
- const MBCartVect pt = 2*xyz - mid;
- return fabs(pt[0]) - dim[0] < etol &&
- fabs(pt[1]) - dim[1] < etol &&
- fabs(pt[2]) - dim[2] < etol &&
- point_in_trilinear_hex( hex, xyz, etol );
-}
-
-
-
-// Wrapper to James Lottes' findpt routines
-// hex_findpt
-// Find the parametric coordinates of an xyz point inside
-// a 3d hex spectral element with n nodes per dimension
-// xm: coordinates fields, value of x,y,z for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
-// n: number of nodes per dimension -- n=2 for a linear element
-// xyz: input, point to find
-// rst: output: parametric coords of xyz inside the element. If xyz is outside the element, rst will be the coords of the closest point
-// dist: output: distance between xyz and the point with parametric coords rst
-extern "C"{
-#include "types.h"
-#include "poly.h"
-#include "tensor.h"
-#include "findpt.h"
-#include "extrafindpt.h"
-#include "errmem.h"
-}
-
-void hex_findpt(real *xm[3],
- int n,
- MBCartVect xyz,
- MBCartVect &rst,
- double &dist)
-{
-
- //compute stuff that only depends on the order -- could be cached
- real *z[3];
- lagrange_data ld[3];
- opt_data_3 data;
-
- //triplicates
- for(int d=0; d<3; d++){
- z[d] = tmalloc(real, n);
- lobatto_nodes(z[d], n);
- lagrange_setup(&ld[d], z[d], n);
- }
-
- opt_alloc_3(&data, ld);
-
- //find nearest point
- real x_star[3];
- xyz.get(x_star);
-
- real r[3] = {0, 0, 0 }; // initial guess for parametric coords
- unsigned c = opt_no_constraints_3;
- dist = opt_findpt_3(&data, (const real **)xm, x_star, r, &c);
- //c tells us if we landed inside the element or exactly on a face, edge, or node
-
- //copy parametric coords back
- rst = r;
-
- //Clean-up (move to destructor if we decide to cache)
- opt_free_3(&data);
- for(int d=0; d<3; ++d)
- lagrange_free(&ld[d]);
- for(int d=0; d<3; ++d)
- free(z[d]);
-}
-
-
-
-
-// hex_eval
-// Evaluate a field in a 3d hex spectral element with n nodes per dimension, at some given parametric coordinates
-// field: field values for each of then n*n*n gauss-lobatto nodes. Nodes are in lexicographical order (x is fastest-changing)
-// n: number of nodes per dimension -- n=2 for a linear element
-// rst: input: parametric coords of the point where we want to evaluate the field
-// value: output: value of field at rst
-
-void hex_eval(real *field,
- int n,
- MBCartVect rstCartVec,
- double &value)
-{
- int d;
- real rst[3];
- rstCartVec.get(rst);
-
- //can cache stuff below
- lagrange_data ld[3];
- real *z[3];
- for(d=0;d<3;++d){
- z[d] = tmalloc(real, n);
- lobatto_nodes(z[d], n);
- lagrange_setup(&ld[d], z[d], n);
- }
-
- //cut and paste -- see findpt.c
- const unsigned
- nf = n*n,
- ne = n,
- nw = 2*n*n + 3*n;
- real *od_work = tmalloc(real, 6*nf + 9*ne + nw);
-
- //piece that we shouldn't want to cache
- for(d=0; d<3; d++){
- lagrange_0(&ld[d], rst[d]);
- }
-
- value = tensor_i3(ld[0].J,ld[0].n,
- ld[1].J,ld[1].n,
- ld[2].J,ld[2].n,
- field,
- od_work);
-
- //all this could be cached
- for(d=0; d<3; d++){
- free(z[d]);
- lagrange_free(&ld[d]);
- }
- free(od_work);
-}
-
-} // namespace MBElemUtil
-
Deleted: MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp 2009-04-17 16:49:56 UTC (rev 2835)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp 2009-04-17 19:44:02 UTC (rev 2836)
@@ -1,38 +0,0 @@
-#include "MBCore.hpp"
-#include "MBCartVect.hpp"
-
-namespace MBElemUtil {
-
- bool nat_coords_trilinear_hex(const MBCartVect* hex_corners,
- const MBCartVect& x,
- MBCartVect& xi,
- double tol);
- bool point_in_trilinear_hex(const MBCartVect *hex_corners,
- const MBCartVect& xyz,
- double etol);
-
- bool point_in_trilinear_hex(const MBCartVect *hex_corners,
- const MBCartVect& xyz,
- const MBCartVect& box_min,
- const MBCartVect& box_max,
- double etol);
-
- //wrapper to hex_findpt
- void nat_coords_trilinear_hex2(const MBCartVect* hex_corners,
- const MBCartVect& x,
- MBCartVect& xi,
- double til);
-
-
-
- void hex_findpt(double *xm[3],
- int n,
- MBCartVect xyz,
- MBCartVect& rst,
- double& dist);
-
- void hex_eval(double *field,
- int n,
- MBCartVect rst,
- double &value);
-} // namespace MBElemUtil
Modified: MOAB/trunk/tools/mbcoupler/Makefile.am
===================================================================
--- MOAB/trunk/tools/mbcoupler/Makefile.am 2009-04-17 16:49:56 UTC (rev 2835)
+++ MOAB/trunk/tools/mbcoupler/Makefile.am 2009-04-17 19:44:02 UTC (rev 2836)
@@ -14,7 +14,6 @@
endif
libmbcoupler_la_SOURCES = \
- MBElemUtil.cpp \
poly.c \
poly.h \
tensor.c \
@@ -29,7 +28,6 @@
MBCoupler.cpp
libmbcoupler_la_include_HEADERS = \
- MBElemUtil.hpp \
MBCoupler.hpp
libmbcoupler_la_includedir = $(includedir)
More information about the moab-dev
mailing list