[MOAB-dev] r2690 - MOAB/trunk/tools/mbcoupler

[acaceres at mcs.anl.gov]

Author: acaceres
Date: 2009-03-10 10:47:23 -0500 (Tue, 10 Mar 2009)
New Revision: 2690

Modified:
   MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
   MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
   MOAB/trunk/tools/mbcoupler/findpt_test.cpp
Log:
Added MBElemUtil:hex_eval, utility function to evaluate the value of a field at given parametric coordinates in a high-order hex element. Wrapper to James Lottes' code.


Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp	2009-03-06 19:08:39 UTC (rev 2689)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp	2009-03-10 15:47:23 UTC (rev 2690)
@@ -180,9 +180,8 @@
 
 
 
-
+// Wrapper to James Lottes' findpt routines
 // hex_findpt
-// Wrapper to James Lottes' findpt routines
 // 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)
@@ -240,4 +239,59 @@
     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
+

Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp	2009-03-06 19:08:39 UTC (rev 2689)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp	2009-03-10 15:47:23 UTC (rev 2690)
@@ -30,4 +30,9 @@
                   MBCartVect xyz, 
                   MBCartVect& rst,
                   double& dist);
+
+  void hex_eval(double *field,
+		int n,
+		MBCartVect rst,
+		double &value);
 } // namespace MBElemUtil

Modified: MOAB/trunk/tools/mbcoupler/findpt_test.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/findpt_test.cpp	2009-03-06 19:08:39 UTC (rev 2689)
+++ MOAB/trunk/tools/mbcoupler/findpt_test.cpp	2009-03-10 15:47:23 UTC (rev 2690)
@@ -51,6 +51,33 @@
 
 
 
+void test_hex_eval()
+{
+    MBCartVect rst(.5,.3,.4);
+    double value;
+
+    const int n=2; //number of nodes per direction (min is 2, for linear element)
+    double *field = tmalloc(double, n*n*n);
+
+    double scale = 1./(n-1);
+    int node = 0;
+
+    //Stuff field with sample data
+    for(int k=0; k<n; k++){
+      for(int j=0; j<n; j++){
+        for(int i=0; i<n; i++){
+          field[node] = 100*scale*(i); 
+          node++;
+        }
+      }
+    }
+        
+    hex_eval(field, n, rst, value);
+    cout << "Value at " << rst << "is: " << value << endl;
+}
+
+
+
 void test_nat_coords_trilinear_hex2()
 {
   MBCartVect hex[8];
@@ -79,7 +106,7 @@
 
 
 int main(){
-
   test_nat_coords_trilinear_hex2();
-    //test_hex_findpt();
+  test_hex_findpt();
+  test_hex_eval();
 }