Templated Compute Kernel Foundation

examples/structural/example_6/example_6.cpp

@@ -51,6 +51,8 @@
 #include "examples/structural/base/inplane_2d_model.h"
 #include "examples/structural/base/truss_2d_model.h"
 #include "examples/structural/base/eyebar_2d_model.h"
+#include "elasticity/elasticity_elem_operations.h"
+#include "base/view.h"
 
 
 // libMesh includes
@@ -1360,7 +1362,7 @@ _optim_con(int*    mode,
 int main(int argc, char* argv[]) {
 
     libMesh::LibMeshInit init(argc, argv);
-
+    
     MAST::Examples::GetPotWrapper
     input(argc, argv, "input");
 
@@ -1401,7 +1403,7 @@ int main(int argc, char* argv[]) {
 
     std::string
     s          = input("optimizer", "optimizer to use in the example", "gcmma");
-
+/*
     if (s == "gcmma") {
 
         optimizer.reset(new MAST::GCMMAOptimizationInterface);
@@ -1447,7 +1449,90 @@ int main(int argc, char* argv[]) {
         else
             optimizer->optimize();
     }
+    */
 
+    TopologyOptimizationSIMP<MAST::Examples::Inplane2DModel>
+    *o = dynamic_cast<TopologyOptimizationSIMP<MAST::Examples::Inplane2DModel>*>(top_opt.get());
+    const libMesh::Elem* e =
+    *(o->_mesh->elements_begin());
+    e->print_info();
+
+    EigenMatrix<Real>::type jac0;
+    {
+        MAST::ElasticityElemOperations<MAST::ElasticityTraits<Real, Real, Real, MAST::ElasticityElemOperationsContext>> elem_ops;
+
+        elem_ops.init(*o->_sys, *o->_discipline, *o->_sys->solution);
+        elem_ops.reinit(*e);
+
+        EigenVector<Real>::type v;
+        EigenMatrix<Real>::type jac;
+        elem_ops.resize_residual(v);
+        elem_ops.resize_jacobian(jac);
+        elem_ops.compute(v, &jac);
+
+        libMesh::out << v << std::endl << jac << std::endl;
+        jac0 = jac;
+    }
+
+    {
+        MAST::ElasticityElemOperations<MAST::ElasticityTraits<Real, Complex, Real, MAST::ElasticityElemOperationsContext>> elem_ops_complex_shape;
+
+        elem_ops_complex_shape.init(*o->_sys, *o->_discipline, *o->_sys->solution);
+        elem_ops_complex_shape.reinit(*e);
+
+        EigenVector<Complex>::type v;
+        EigenMatrix<Complex>::type jac;
+        EigenMatrix<Real>::type jac1;
+        elem_ops_complex_shape.resize_residual(v);
+        elem_ops_complex_shape.resize_jacobian(jac);
+        elem_ops_complex_shape.compute(v, &jac);
+        elem_ops_complex_shape.compute_complex_step_jacobian<Real>(jac1);
+
+        libMesh::out << v << std::endl << jac << std::endl << jac1 << std::endl;
+        jac1 -= jac0;
+        libMesh::out << jac1 << std::endl << jac1.norm() << std::endl;
+    }
+
+    {
+        MAST::ElasticityElemOperations<MAST::ElasticityTraits<Real, Real, Complex, MAST::ElasticityElemOperationsContext>> elem_ops_complex_sol;
+
+        elem_ops_complex_sol.init(*o->_sys, *o->_discipline, *o->_sys->solution);
+        elem_ops_complex_sol.reinit(*e);
+
+        EigenVector<Complex>::type v;
+        EigenMatrix<Complex>::type jac;
+        EigenMatrix<Real>::type jac1;
+
+        elem_ops_complex_sol.resize_residual(v);
+        elem_ops_complex_sol.resize_jacobian(jac);
+        elem_ops_complex_sol.compute(v, &jac);
+        elem_ops_complex_sol.compute_complex_step_jacobian<Real>(jac1);
+
+        libMesh::out << v << std::endl << jac << std::endl << jac1 << std::endl;
+        jac1 -= jac0;
+        libMesh::out << jac1 << std::endl << jac1.norm() << std::endl;
+    }
+
+    {
+        MAST::ElasticityElemOperations<MAST::ElasticityTraits<Real, Complex, Complex, MAST::ElasticityElemOperationsContext>> elem_ops_complex_sol2;
+
+        elem_ops_complex_sol2.init(*o->_sys, *o->_discipline, *o->_sys->solution);
+        elem_ops_complex_sol2.reinit(*e);
+
+        EigenVector<Complex>::type v;
+        EigenMatrix<Complex>::type jac;
+        EigenMatrix<Real>::type jac1;
+
+        elem_ops_complex_sol2.resize_residual(v);
+        elem_ops_complex_sol2.resize_jacobian(jac);
+        elem_ops_complex_sol2.compute(v, &jac);
+        elem_ops_complex_sol2.compute_complex_step_jacobian<Real>(jac1);
+
+        libMesh::out << v << std::endl << jac << std::endl << jac1 << std::endl;
+        jac1 -= jac0;
+        libMesh::out << jac1 << std::endl << jac1.norm() << std::endl;
+    }
+
     // END_TRANSLATE
     return 0;
 }

src/base/CMakeLists.txt

@@ -11,6 +11,7 @@ target_sources(mast
         ${CMAKE_CURRENT_LIST_DIR}/complex_assembly_elem_operations.h
         ${CMAKE_CURRENT_LIST_DIR}/complex_mesh_field_function.cpp
         ${CMAKE_CURRENT_LIST_DIR}/complex_mesh_field_function.h
+        ${CMAKE_CURRENT_LIST_DIR}/compute_kernel.h
         ${CMAKE_CURRENT_LIST_DIR}/constant_field_function.cpp
         ${CMAKE_CURRENT_LIST_DIR}/constant_field_function.h
         ${CMAKE_CURRENT_LIST_DIR}/dof_coupling_base.cpp

src/base/computation_base.h

@@ -0,0 +1,34 @@
+
+#ifndef _mast_computation_base_h_
+#define _mast_computation_base_h_
+
+// MAST includes
+#include "base/compute_kernel_base.h"
+
+namespace MAST {
+
+/*! Collection of compute kernels with specified dependencies and data views*/
+template <typename Traits, typename ContextType>
+class ComputationBase {
+
+public:
+
+    using basis_scalar_type = typename Traits::basis_scalar_type;
+    using nodal_scalar_type = typename Traits::nodal_scalar_type;
+    using sol_scalar_type   = typename Traits::sol_scalar_type;
+
+    ComputationBase() { }
+    virtual ~ComputationBase() { }
+    void add_compute_kernel(MAST::ComputeKernelBase<ContextType>& c);
+    /*! parses through all the compute kernels and their views to prepare a graph of dependency. */
+    void prepare();
+    void print_graph();
+    void execute();
+
+protected:
+
+};
+
+}
+
+#endif // _mast_computation_base_h_

src/base/compute_kernel.h

@@ -0,0 +1,161 @@
+//
+//  compute_kernel.h
+//  compute_kernels
+//
+//  Created by Manav Bhatia on 5/1/20.
+//  Copyright © 2020 Manav Bhatia. All rights reserved.
+//
+
+#ifndef _mast_compute_kernel_h_
+#define _mast_compute_kernel_h_
+
+// C++ includes
+#include <vector>
+#include <string>
+
+// MAST includes
+#include "base/compute_kernel_base.h"
+#include "base/view.h"
+
+
+namespace MAST {
+
+
+template <typename ValueType>
+struct IsScalarType {
+    using type = std::false_type;
+};
+
+
+template <>
+struct IsScalarType<Real> {
+    using type = std::true_type;
+};
+
+
+template <>
+struct IsScalarType<std::complex<Real>> {
+    using type = std::true_type;
+};
+
+
+
+template <typename KernelValueType, typename IsIndexable, typename IsScalarType = typename MAST::IsScalarType<KernelValueType>::type>
+struct KernelReturnType {
+
+};
+
+
+template <typename KernelValueType>
+struct KernelReturnType<KernelValueType, std::false_type, std::false_type> {
+
+    using type = KernelValueType;
+};
+
+
+template <typename KernelValueType>
+struct KernelReturnType<KernelValueType, std::false_type, std::true_type> {
+
+    using type = KernelValueType;
+};
+
+
+
+template <typename KernelValueType>
+struct KernelReturnType<KernelValueType,  std::true_type, std::false_type> {
+
+    using type = Eigen::Map<KernelValueType>;
+};
+
+
+template <typename KernelValueType>
+struct KernelReturnType<KernelValueType,  std::true_type, std::true_type> {
+
+    using type = KernelValueType&;
+};
+
+
+template <typename KernelType, typename KernelViewType = typename KernelType::view_type, typename ContextType>
+KernelViewType build_kernel_view(KernelType& k, ContextType& c) { }
+
+
+
+
+// Forward decleration
+template <typename> class ComputeKernelDerivative;
+
+
+template <typename ContextType>
+class ComputeKernel: public MAST::ComputeKernelBase<ContextType> {
+
+public:
+
+    using derivative_kernel_type = MAST::ComputeKernelDerivative<ContextType>;
+
+    ComputeKernel(const std::string& nm, const bool executable):
+    MAST::ComputeKernelBase<ContextType>(nm, executable),
+    _derivative_kernel   (nullptr)
+    {}
+
+    virtual ~ComputeKernel() {}
+
+    virtual inline void set_derivative_kernel(derivative_kernel_type& d)
+    {
+        libmesh_assert_msg(!_derivative_kernel, "Derivative kernel already set.");
+        _derivative_kernel = &d;
+    }
+
+    virtual inline const derivative_kernel_type& get_derivative_kernel() const
+    {
+        libmesh_assert_msg(_derivative_kernel, "Derivative kernel not set");
+        return *_derivative_kernel;
+    }
+
+protected:
+
+    derivative_kernel_type     *_derivative_kernel;
+};
+
+
+
+template <typename ContextType>
+class ComputeKernelDerivative: public MAST::ComputeKernelBase<ContextType> {
+
+public:
+
+    using primal_kernel_type   = MAST::ComputeKernel<ContextType>;
+
+    ComputeKernelDerivative (const std::string& nm, const bool executable):
+    MAST::ComputeKernelBase<ContextType> (nm),
+    _primal_kernel          (nullptr),
+    _f                      (nullptr)
+    {}
+
+    virtual ~ComputeKernelDerivative() {}
+
+    virtual inline void  set_primal_kernel(primal_kernel_type& k)
+    {
+        libmesh_assert_msg(!_primal_kernel, "Primal kernel already set.");
+        _primal_kernel = &k;
+    }
+
+    virtual inline primal_kernel_type& get_primal_kernel()
+    {
+        libmesh_assert_msg(_primal_kernel, "Primal kernel not set.");
+        return *_primal_kernel;
+    }
+
+    virtual inline void  set_derivative_paramter(const MAST::FunctionBase& f);
+
+protected:
+
+    primal_kernel_type *_primal_kernel;
+    const MAST::FunctionBase       *_f;
+};
+
+
+
+}
+
+#endif /* __mast_compute_kernel_h__ */
+

src/base/compute_kernel_base.h

@@ -0,0 +1,59 @@
+
+#ifndef _mast_compute_kernel_base_h_
+#define _mast_compute_kernel_base_h_
+
+// C++ includes
+#include <set>
+#include <memory>
+
+// MAST includes
+#include "base/mast_data_types.h"
+
+
+namespace MAST {
+
+// Forward declerations
+class FunctionBase;
+
+template <typename ContextType>
+class ComputeKernelBase {
+
+public:
+
+    ComputeKernelBase(const std::string& nm,
+                      const bool executable):
+    _nm          (nm),
+    _executable  (executable) {}
+
+    virtual ~ComputeKernelBase() {}
+
+    virtual inline bool is_executable() const { return _executable;}
+    virtual inline bool depends_on(const MAST::ComputeKernelBase<ContextType>& d) const
+    { return _dependency.count(&d);}
+    virtual inline const std::set<const MAST::ComputeKernelBase<ContextType>*>& get_dependencies() const
+    { return _dependency;}
+
+protected:
+
+    virtual inline void _add_dependency(const MAST::ComputeKernelBase<ContextType>& d) { _dependency.insert(&d);}
+
+    const std::string _nm;
+    const bool        _executable;
+    std::set<const MAST::ComputeKernelBase<ContextType>*> _dependency;
+};
+
+
+struct CurrentComputation {
+
+    CurrentComputation(): time (0.), dt(0.), param(nullptr) {}
+
+    std::vector<uint_type>  current_indices;
+    Real                    time;
+    Real                    dt;
+    /*! parameter for which sensitivity is being computed */
+    FunctionBase     *param;
+};
+
+}
+
+#endif // _mast_compute_kernel_base_h_