Feat/add common finite element utilities

CMakeLists.txt

@@ -16,7 +16,7 @@ set(TARDIGRADE_BALANCE_EQUATIONS_USE_LIBXSMM OFF CACHE BOOL "Flag for whether to
 
 # Set the internal support libraries
 set(INTERNAL_SUPPORT_LIBRARIES )
-set(ADDITIONAL_HEADER_ONLY_LIBRARIES "tardigrade_balance_of_mass")
+set(ADDITIONAL_HEADER_ONLY_LIBRARIES "tardigrade_balance_of_mass" "tardigrade_finite_element_utilities")
 set(PROJECT_SOURCE_FILES ${PROJECT_NAME}.cpp ${PROJECT_NAME}.h)
 set(PROJECT_PRIVATE_HEADERS "")
 foreach(package ${INTERNAL_SUPPORT_LIBRARIES})

docs/sphinx/changelog.rst

@@ -14,3 +14,4 @@ New Features
 - Initial commit of the balance equation repository (:pull:`1`). By `Nathan Miller`_.
 - Removed mentions of tardigrade hydra from the readme (:pull:`2`). By `Nathan Miller`_.
 - Added the calculation of the balance of mass (:pull:`3`). By `Nathan Miller`_.
+- Added the calculation of the derivative of the spatial gradient of a quantity (:pull:`4`). By `Nathan Miller`_.

src/cpp/tardigrade_balance_of_mass.h

@@ -6,10 +6,9 @@
   ******************************************************************************
   */
 
-#ifndef TARDIGRADE_MASS_CHANGE_DEFORMATION_H
-#define TARDIGRADE_MASS_CHANGE_DEFORMATION_H
+#ifndef TARDIGRADE_BALANCE_OF_MASS_H
+#define TARDIGRADE_BALANCE_OF_MASS_H
 
-#include<vector>
 #include<array>
 
 #include "tardigrade_error_tools.h"

src/cpp/tardigrade_finite_element_utilities.cpp

@@ -0,0 +1,53 @@
+/**
+  ******************************************************************************
+  * \file tardigrade_finite_element_utilities.cpp
+  ******************************************************************************
+  * The source file for utilities which can assist using the balance equations
+  * in finite element codes. We here assume that unknown quantities are
+  * computed at the evaluation point using interpolation functions (here called
+  * interp) and are projected to the nodes using test functions (here called
+  * test).
+  ******************************************************************************
+  */
+
+namespace tardigradeBalanceEquations{
+
+    namespace finiteElementUtilities{
+
+        template<typename grad_iterator, typename output_iterator>
+        void computeGradientSpatialJacobian( const grad_iterator &grad_a_start, const unsigned int grad_a_size,
+                                             floatVector grad_test, const unsigned int index, output_iterator dgrad_adui_start ){
+            /*!
+             * Compute the derivative of the spatial gradient of a quantity a
+             * in the current configuration w.r.t. the spatial degrees of freedom
+             * 
+             * \f$ \frac{D}{Du_a} \left( a_{i,j} \right) \f$
+             * 
+             * \param &grad_a_start: An iterator representing the start of the gradient of the quantity
+             * \param &grad_a_size: The size of the gradient of a
+             * \param &grad_test: The gradient of the test function
+             * \param &index: The index of the spatial degree of freedom to compute the derivative for (0, 1 or 2)
+             * \param &dgrad_adui_start: An iterator representing the start of the output
+             */
+
+            TARDIGRADE_ERROR_TOOLS_CHECK( ( grad_a_size % dim ) == 0, "The incoming spatial gradient has a dimension of " + std::to_string( grad_a_size ) + " which is not a multiple of " + std::to_string( dim ) );
+
+            const unsigned int a_size = grad_a_size / dim;
+
+            std::fill( dgrad_adui_start, dgrad_adui_start + a_size, 0 );
+
+            for ( unsigned int i = 0; i < a_size; i++ ){
+
+                for ( unsigned int j = 0; j < dim; j++ ){
+
+                    *( dgrad_adui_start + dim * i + j ) -= *( grad_a_start + dim * i + index ) * grad_test[ j ];
+
+                }
+
+            }
+
+        }
+
+    }
+
+}

src/cpp/tardigrade_finite_element_utilities.h

@@ -0,0 +1,44 @@
+/**
+  ******************************************************************************
+  * \file tardigrade_finite_element_utilities.h
+  ******************************************************************************
+  * The header file for utilities which can assist using the balance equations
+  * in finite element codes. We here assume that unknown quantities are
+  * computed at the evaluation point using interpolation functions (here called
+  * interp) and are projected to the nodes using test functions (here called
+  * test).
+  ******************************************************************************
+  */
+
+#ifndef TARDIGRADE_FINITE_ELEMENT_UTILITIES_H
+#define TARDIGRADE_FINITE_ELEMENT_UTILITIES_H
+
+#include<array>
+
+#include "tardigrade_error_tools.h"
+
+namespace tardigradeBalanceEquations{
+
+    namespace finiteElementUtilities{
+
+        constexpr unsigned int dim = 3; //!< Set the dimension as 3D by default
+
+        constexpr unsigned int sot_dim = dim * dim; //!< Set the dimensions of a standard second order tensor
+
+        typedef double floatType; //!< Define the float type as a double
+
+        typedef std::array< floatType, dim > floatVector; //!< Define a standard vector
+
+        typedef std::array< floatType, sot_dim > secondOrderTensor; //!< Define a standard second-order tensor
+
+        template<typename grad_iterator, typename output_iterator>
+        void computeGradientSpatialJacobian( const grad_iterator &grad_a_start, const unsigned int grad_a_size, 
+                                             floatVector grad_test, const unsigned int index, output_iterator dgrad_adui_start );
+
+    }
+
+}
+
+#include "tardigrade_finite_element_utilities.cpp"
+
+#endif

src/cpp/tests/test_tardigrade_finite_element_utilities.cpp

@@ -0,0 +1,77 @@
+/**
+  * \file test_tardigrade_finite_element_utilities.cpp
+  *
+  * Tests for tardigrade_finite_element_utilities
+  */
+
+#include<tardigrade_finite_element_utilities.h>
+#include<sstream>
+#include<fstream>
+#include<iostream>
+#include<array>
+
+#define BOOST_TEST_MODULE test_tardigrade_finite_element_utilities
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/tools/output_test_stream.hpp>
+
+#define DEFAULT_TEST_TOLERANCE 1e-6
+#define CHECK_PER_ELEMENT boost::test_tools::per_element( )
+
+struct cout_redirect{
+    cout_redirect( std::streambuf * new_buffer )
+        : old( std::cout.rdbuf( new_buffer ) )
+    { }
+
+    ~cout_redirect( ) {
+        std::cout.rdbuf( old );
+    }
+
+    private:
+        std::streambuf * old;
+};
+
+struct cerr_redirect{
+    cerr_redirect( std::streambuf * new_buffer )
+        : old( std::cerr.rdbuf( new_buffer ) )
+    { }
+
+    ~cerr_redirect( ) {
+        std::cerr.rdbuf( old );
+    }
+
+    private:
+        std::streambuf * old;
+};
+
+typedef tardigradeBalanceEquations::finiteElementUtilities::floatType floatType; //!< Define the float type to be the same as in the balance of mass
+
+typedef tardigradeBalanceEquations::finiteElementUtilities::floatVector floatVector; //!< Define the float vector type to be the same as in the balance of mass
+
+typedef tardigradeBalanceEquations::finiteElementUtilities::secondOrderTensor secondOrderTensor; //!< Define the second order tensor type to be the same as in the balance of mass
+
+BOOST_AUTO_TEST_CASE( test_computeGradientSpatialJacobian, * boost::unit_test::tolerance( DEFAULT_TEST_TOLERANCE ) ){
+    /*!
+     * Test the computation of the jacobian of the gradient of a quantity w.r.t. the spatial displacement
+     */
+
+    std::array< floatType, 12 > grad_a = { 0,  1,  2,
+                                           3,  4,  5,
+                                           6,  7,  8,
+                                           9, 10, 11 };
+
+    floatVector grad_test = { 0.1, 0.2, 0.3 };
+
+    std::array< floatType, 12 > answer = { -0.2, -0.4, -0.6,
+                                           -0.5, -1. , -1.5,
+                                           -0.8, -1.6, -2.4,
+                                           -1.1, -2.2, -3.3 };
+
+    std::array< floatType, 12 > result;
+    std::fill( std::begin( result ), std::end( result ), 0 );
+
+    tardigradeBalanceEquations::finiteElementUtilities::computeGradientSpatialJacobian( std::begin( grad_a ), 12, grad_test, 2, std::begin( result ) );
+
+    BOOST_TEST( answer == result, CHECK_PER_ELEMENT );
+
+}
+