Merging develop into master

Author: Nightly Bot

CMakeLists.txt

@@ -136,6 +136,7 @@ add_subdirectory(pumi)
 add_subdirectory(ma)
 add_subdirectory(crv)
 add_subdirectory(spr)
+add_subdirectory(ree)
 add_subdirectory(sam)
 add_subdirectory(phasta)
 add_subdirectory(stk)

README.md

@@ -26,6 +26,7 @@ For more information, start at our
 * MTH: Math containers and routines
 * CRV: Support for curved meshes with Bezier Shapes
 * PYCORE: Python Wrappers (see python_wrappers/README.md for build instructions)
+* REE: Residual based implicit error estimator
 
 ### How do I get set up? ###
 

apf/apf.cc

@@ -4,7 +4,6 @@
  * This work is open source software, licensed under the terms of the
  * BSD license as described in the LICENSE file in the top-level directory.
  */
-
 #include "apf.h"
 #include "apfScalarField.h"
 #include "apfScalarElement.h"
@@ -464,28 +463,20 @@ void getVectorShapeValues(Element* e, Vector3 const& local,
   NewArray<Vector3> vvals(values.size());
   e->getShape()->getVectorValues(e->getMesh(), e->getEntity(), local, vvals);
 
-  // Perform Piola transformation
-  if( e->getDimension() == e->getMesh()->getDimension() ) // i.e. J is square
-  {
-    apf::Matrix3x3 Jinv;
-    apf::getJacobianInv( e->getParent(), local, Jinv );
-    apf::Matrix3x3 JinvT = apf::transpose(Jinv);
-
-    // u(x_hat) * J(x_hat)^{-1}
-    for( size_t i = 0; i < values.size(); i++ ) {
-      for ( int j = 0; j < 3; j++ ) {
-        values[i][j] = 0.;
-        for ( int k = 0; k < 3; k++ )
-          values[i][j] += vvals[i][k] * JinvT[k][j];
-      }
+  apf::Matrix3x3 Jinv;
+  apf::getJacobianInv( e->getParent(), local, Jinv );
+  apf::Matrix3x3 JinvT = apf::transpose(Jinv);
+
+  // Perform Piola transformation - u(x_hat) * J(x_hat)^{-1}
+  int d = 0;
+  (e->getDimension() == e->getMesh()->getDimension()) ? d = 3 : d = 2;
+  for( size_t i = 0; i < values.size(); i++ ) {
+    for ( int j = 0; j < 3; j++ ) {
+      values[i][j] = 0.;
+      for ( int k = 0; k < d; k++ )
+        values[i][j] += vvals[i][k] * JinvT[k][j];
     }
   }
-  else
-  {
-    // TODO when ref dim != mesh space dim. Pseudo-inverse needed.
-    PCU_ALWAYS_ASSERT_VERBOSE(false,
-    	"not yet implemented for 3D surface meshes (i.e., manifolds)!");
-  }
 }
 
 void getCurlShapeValues(Element* e, Vector3 const& local,

apf/apfElement.cc

@@ -119,4 +119,11 @@ void Element::getNodeData()
   field->getData()->getElementData(entity,nodeData);
 }
 
+void Element::getElementNodeData(NewArray<double>& d)
+{
+  d.allocated() ? d.resize(nen) : d.allocate(nen);
+  for (int i = 0; i < nen; i++)
+    d[i] = nodeData[i];
+}
+
 }//namespace apf

apf/apfElement.h

@@ -35,6 +35,7 @@ class Element
     EntityShape* getShape() {return shape;}
     FieldShape* getFieldShape() {return field->getShape();}
     void getComponents(Vector3 const& xi, double* c);
+    void getElementNodeData(NewArray<double>& d);
   protected:
     void init(Field* f, MeshEntity* e, VectorElement* p);
     void getNodeData();

mth/mth_def.h

@@ -57,12 +57,12 @@ template <class T, unsigned M, unsigned N>
 void transpose(Matrix<T,M,N> const& a,
     Matrix<T,N,M>& b)
 {
-  unsigned m = a.rows();
+  unsigned m = a.cols();
   unsigned n = a.rows();
   b.resize(m, n);
   for (unsigned i=0; i < m; ++i)
-  for (unsigned j=0; j < n; ++j)
-    b(j,i) = a(i,j);
+    for (unsigned j=0; j < n; ++j)
+      b(i,j) = a(j,i);
 }
 
 template <class T>

pumi-meshes

@@ -0,0 +1,32 @@
+if(DEFINED TRIBITS_PACKAGE)
+  include(pkg_tribits.cmake)
+  return()
+endif()
+
+# Package sources
+set(SOURCES
+  reeResidualFunctionals.cc
+  reeFluxCorrection.cc
+  reeCorrectedFlux.cc
+  reeEstimateError.cc
+  reeSizeField.cc)
+
+# Package headers
+set(HEADERS
+  ree.h)
+
+# Add the ree library
+add_library(ree ${SOURCES})
+
+# Include directories
+target_include_directories(ree INTERFACE
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+    $<INSTALL_INTERFACE:include>
+    )
+
+# Link this package to these libraries
+target_link_libraries(ree PUBLIC apf pcu)
+
+scorec_export_library(ree)
+
+bob_end_subdir()

ree/CMakeLists.txt

@@ -0,0 +1,3 @@
+TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
+  LIB_REQUIRED_PACKAGES SCORECapf
+  )

ree/cmake/Dependencies.cmake

@@ -0,0 +1,23 @@
+tribits_package(SCORECree)
+
+# THIS IS WHERE TRIBITS GETS HEADERS
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+#Sources & Headers
+set(SOURCES
+  reeResidualFunctionals.cc
+  reeFluxCorrection.cc
+  reeCorrectedFlux.cc
+  reeEstimateError.cc
+  reeSizeField.cc)
+
+set(HEADERS
+  ree.h)
+
+#Library
+tribits_add_library(
+  ree
+  HEADERS ${HEADERS}
+  SOURCES ${SOURCES})
+
+tribits_package_postprocess()