From a440335bc53cc802fd2ec580d30471856d12a78c Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Wed, 13 May 2020 00:49:06 -0500 Subject: [PATCH 1/9] -- First commit of a templated compute kernel foundation to enable complex-step and auto-diff sensitivity analysis. -- Skeleton implementation of strain energy kernel for elasticity computations and a code to set this up. These classes are instantiated for different combinations of variable types in example 6. --- examples/structural/example_6/example_6.cpp | 6 + src/base/compute_kernel.h | 804 ++++++++++++++++++++ 2 files changed, 810 insertions(+) create mode 100644 src/base/compute_kernel.h diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp index 3d885464..14b1c3d0 100644 --- a/examples/structural/example_6/example_6.cpp +++ b/examples/structural/example_6/example_6.cpp @@ -51,6 +51,7 @@ #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 "base/compute_kernel.h" // libMesh includes @@ -1448,6 +1449,11 @@ int main(int argc, char* argv[]) { optimizer->optimize(); } + MAST::ElasticityElemOperations> elem_ops; + MAST::ElasticityElemOperations> elem_ops_complex_shape; + MAST::ElasticityElemOperations> elem_ops_complex_sol; + MAST::ElasticityElemOperations> elem_ops_complex_sol2; + // END_TRANSLATE return 0; } diff --git a/src/base/compute_kernel.h b/src/base/compute_kernel.h new file mode 100644 index 00000000..73405888 --- /dev/null +++ b/src/base/compute_kernel.h @@ -0,0 +1,804 @@ +// +// 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 +#include + +// MAST includes +#include "base/mast_data_types.h" + +// libMesh includes +#include "libmesh/quadrature.h" +#include "libmesh/fe_base.h" +#include "libmesh/elem.h" + +typedef unsigned int uint_type; + +namespace MAST { + + +// Forward declerations +class FunctionBase; + + +class ComputeKernelBase { + +public: + + ComputeKernelBase(const std::string& nm): _nm(nm) {} + virtual ~ComputeKernelBase() {} + virtual inline bool depends_on(const MAST::ComputeKernelBase& d) const { return _dependency.count(&d);} + virtual inline const std::set& get_dependencies() const { return _dependency;} + virtual inline void pre_execute() {} + virtual inline void post_execute() {} + virtual inline void execute() = 0; + +protected: + + virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} + + const std::string _nm; + std::set _dependency; +}; + + +struct CurrentComputation { + + CurrentComputation(): qp(0), time (0.), dt(0.), param(nullptr) {} + + uint_type qp; + Real time; + Real dt; + /*! parameter for which sensitivity is being computed */ + FunctionBase *param; +}; + + +template +class Quadrature: public MAST::ComputeKernelBase { + +public: + + using scalar_type = ScalarType; + + Quadrature(const std::string& nm): MAST::ComputeKernelBase(nm) {} + virtual ~Quadrature() {} + virtual inline uint_type dim() const = 0; + virtual inline uint_type order() const = 0; + virtual inline uint_type n_points() const = 0; + virtual inline ScalarType qp_coord(uint_type qp, uint_type xi_i) const = 0; + virtual inline ScalarType weight(uint_type qp) const = 0; + +}; + + + +/*! This provides a */ +template +class MappedQuadrature: public MAST::Quadrature { + +public: + + using scalar_type = ScalarType; + + MappedQuadrature(const std::string& nm): MAST::Quadrature(nm), _dim(0) {} + virtual ~MappedQuadrature() {} + virtual void set_data(const uint_type dim, + const typename ViewTraits::quadrature_point_type points, + const typename ViewTraits::quadrature_point_type weights) { + _points = points; + _weights = weights; + } + virtual inline uint_type dim() const override { return _dim;} + virtual inline uint_type order() const override { /* should not be called*/ libmesh_assert(false);} + virtual inline uint_type n_points() const override { return _weights.size();} + virtual inline ScalarType qp_coord(uint_type qp, uint_type xi_i) const override { return _points(qp, xi_i);} + virtual inline ScalarType weight(uint_type qp) const { return _weights(qp);} + +protected: + + uint_type _dim; + typename ViewTraits::quadrature_point_type _points; + typename ViewTraits::quadrature_point_type _weights; +}; + + + +/*! serves as a wrapper around libMesh */ +class libMeshQuadrature: public MAST::Quadrature { + +public: + + /*! the quadrature object \p q is expected to be initialized outside of this class. */ + libMeshQuadrature(const std::string& nm): + MAST::Quadrature(nm), + _q (nullptr) + {} + virtual ~libMeshQuadrature() {} + virtual inline void execute() override { } + const libMesh::QBase& get_libmesh_object() const { return *_q;} + virtual inline uint_type dim() const override { return _q->get_dim();} + virtual inline uint_type order() const override { return _q->get_order();} + virtual inline uint_type n_points() const override { return _q->n_points();} + virtual inline scalar_type qp_coord(uint_type qp, uint_type xi_i) const override + { return _q->get_points()[qp](xi_i);} + virtual inline scalar_type weight(uint_type qp) const override { return _q->w(qp);} + +protected: + + libMesh::QBase* _q; +}; + + + +template +class FEBasis: ComputeKernelBase { + +public: + + using basis_scalar_type = ScalarType; + + FEBasis(const std::string& nm): + MAST::ComputeKernelBase(nm), + _quadrature (nullptr) + { } + virtual ~FEBasis() {} + + virtual inline void set_quadrature(const MAST::Quadrature& q) { _quadrature = &q;} + virtual inline uint_type dim() const = 0; + virtual inline uint_type order() const = 0; + virtual inline uint_type n_basis() const = 0; + virtual inline ScalarType qp_coord(uint_type qp, uint_type x_i) const + { return _quadrature->qp_coord(qp, x_i); } + virtual inline ScalarType phi(uint_type qp, uint_type phi_i) const = 0; + virtual inline ScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + +protected: + + const MAST::Quadrature *_quadrature; +}; + + +/*! serves as a wrapper around libMesh FEBase object to provide */ +class libMeshFE: public FEBasis { + +public: + + using scalar_type = Real; + + libMeshFE(const std::string& nm): + MAST::FEBasis(nm), + _compute_dphi_dxi (false), + _fe (nullptr) + { } + virtual ~libMeshFE() { } + + inline void set_compute_dphi_dxi(bool f) { _compute_dphi_dxi = f;} + + virtual inline void reinit(const libMesh::FEBase& fe) { + _fe = &fe; + const libMesh::FEMap& m = fe.get_fe_map(); + _dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; + } + + virtual inline void execute() override { } + + virtual inline uint_type dim() const override { return _fe->get_dim();} + + virtual inline uint_type order() const override {return _fe->get_order();} + + virtual inline uint_type n_basis() const override { return _fe->n_shape_functions();} + + virtual inline basis_scalar_type phi(uint_type qp, uint_type phi_i) const override + {return _fe->get_phi()[phi_i][qp];} + + virtual inline basis_scalar_type + dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const override + { return (*_dphi_dxi[xi_i])[phi_i][qp];} + +protected: + + bool _compute_dphi_dxi; + const libMesh::FEBase* _fe; + std::vector>*> _dphi_dxi; +}; + + +template +struct DeducedScalarType { }; + +template <> +struct DeducedScalarType { using type = Real;}; + +template <> +struct DeducedScalarType, Real> { using type = std::complex;}; + +template <> +struct DeducedScalarType> { using type = std::complex;}; + +template <> +struct DeducedScalarType, std::complex> { using type = std::complex;}; + + + +template +struct EigenFEShapeDataViewTraits { + + using xyz_view_type = Eigen::Matrix; + using detJ_view_type = Eigen::Matrix; + using detJxW_view_type = Eigen::Matrix; + using dx_dxi_view_type = Eigen::Matrix; + using dxi_dx_view_type = Eigen::Matrix; + using dphi_dx_view_type = Eigen::Matrix; + + template + void init_view(ViewType& view, uint_type n1); + + template + void init_view(ViewType& view, uint_type n1, uint_type n2); +}; + + + +template +struct EigenFESolDataViewTraits { + + using coefficient_view_type = Eigen::Matrix; + using u_view_type = Eigen::Matrix; + using du_dx_view_type = Eigen::Matrix; + + template + void init_view(ViewType& view, uint_type n1); + + template + void init_view(ViewType& view, uint_type n1, uint_type n2); +}; + + + +template +class FEShapeDataBase: public MAST::ComputeKernelBase { + +public: + + FEShapeDataBase(const std::string& nm): + MAST::ComputeKernelBase(nm), + _compute_xyz (false), + _compute_Jac (false), + _compute_detJ (false), + _compute_JxW (false), + _compute_dphi_dx (false), + _compute_normal (false), + _spatial_dim (0), + _fe_basis (nullptr) + {} + virtual ~FEShapeDataBase() {} + inline void set_compute_xyz(bool f) { _compute_xyz = f;} + inline void set_compute_Jac(bool f) { _compute_Jac = f;} + inline void set_compute_detJ(bool f) { _compute_detJ = f;} + inline void set_compute_detJxW(bool f) { _compute_JxW = f;} + inline void set_compute_dphi_dx(bool f) { _compute_dphi_dx = f;} + inline void set_compute_normal(bool f) { _compute_normal = f;} + inline void set_fe_basis(MAST::FEBasis& basis) { _fe_basis = &basis;} + virtual inline void execute() override { } + //virtual inline void reinit_for_side(MAST::FEBasis& basis, uint_type s) = 0; + virtual inline uint_type ref_dim() const { return _fe_basis->dim();} + virtual inline uint_type spatial_dim() const { return _spatial_dim;} + virtual inline uint_type order() const { return _fe_basis->order();} + virtual inline uint_type n_basis() const { return _fe_basis->n_basis();} + virtual inline BasisScalarType phi(uint_type qp, uint_type phi_i) const + { return _fe_basis->phi(qp, phi_i);} + virtual inline BasisScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const + { return _fe_basis->dphi_dxi(qp, phi_i, xi_i);} + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const = 0; + virtual inline NodalScalarType detJ(uint_type qp) const = 0; + virtual inline NodalScalarType detJxW(uint_type qp) const = 0; + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + +protected: + + bool _compute_xyz; + bool _compute_Jac; + bool _compute_detJ; + bool _compute_JxW; + bool _compute_dphi_dx; + bool _compute_normal; + uint_type _spatial_dim; + + MAST::FEBasis *_fe_basis; +}; + + + +/*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate + nodal locations. Typically, Lagrange shape functoins are used for this purpose */ +template +class FEGeometryBasisDerivedData: +public MAST::FEShapeDataBase { + +public: + + FEGeometryBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm) {} + virtual ~FEGeometryBasisDerivedData() {} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dx_dxi(qp, x_i, xi_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dxi_dx(qp, x_i, xi_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, phi_i, x_i);} + +protected: + + typename ViewTraits::xyz_view_type _xyz; + typename ViewTraits::detJ_view_type _detJ; + typename ViewTraits::detJxW_view_type _detJxW; + typename ViewTraits::dx_dxi_view_type _dx_dxi; + typename ViewTraits::dxi_dx_view_type _dxi_dx; + typename ViewTraits::dphi_dx_view_type _dphi_dx; +}; + + + +/*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate + nodal locations. Typically, Lagrange shape functoins are used for this purpose */ +template +class FEGeometryBasisDerivedData>: +public MAST::FEShapeDataBase { + +public: + + FEGeometryBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm) {} + + virtual ~FEGeometryBasisDerivedData() {} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} + +protected: + + typename EigenFEShapeDataViewTraits::xyz_view_type _xyz; + typename EigenFEShapeDataViewTraits::detJ_view_type _detJ; + typename EigenFEShapeDataViewTraits::detJxW_view_type _detJxW; + typename EigenFEShapeDataViewTraits::dx_dxi_view_type _dx_dxi; + typename EigenFEShapeDataViewTraits::dxi_dx_view_type _dxi_dx; + typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; +}; + + + +/*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, + this is any non-Lagrange basis.*/ +template +class FEBasisDerivedData: public MAST::FEShapeDataBase { + +public: + + using fe_geom_derived_type = MAST::FEGeometryBasisDerivedData; + + FEBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm), + _fe_geom(nullptr) + {} + + virtual ~FEBasisDerivedData() {} + + virtual void inline set_fe_geom_derived_data(const fe_geom_derived_type& fe_geom) + { _fe_geom = &fe_geom;} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _fe_geom->xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _fe_geom->detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _fe_geom->detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dx_dxi(qp, x_i, xi_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dxi_dx(qp, x_i, xi_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, phi_i, x_i);} + +protected: + + const fe_geom_derived_type *_fe_geom; + + typename ViewTraits::dphi_dx_view_type _dphi_dx; +}; + + + +/*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, + this is any non-Lagrange basis.*/ +template +class FEBasisDerivedData>: +public MAST::FEShapeDataBase { + +public: + + using fe_derived_data_type = MAST::FEGeometryBasisDerivedData>; + + FEBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm), + _fe_geom (nullptr) + {} + + virtual ~FEBasisDerivedData() {} + + virtual void inline set_fe_geom_derived_data(const fe_derived_data_type& fe_geom) + { _fe_geom = &fe_geom;} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _fe_geom->xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _fe_geom->detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _fe_geom->detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dx_dxi(qp, x_i, xi_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dxi_dx(qp, x_i, xi_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} + +protected: + + const fe_derived_data_type& _fe_geom; + + typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; +}; + + + +/*! provides access to the element solution vector through a memory view. */ +template +class FieldVariable { + +public: + + FieldVariable(); + virtual ~FieldVariable(); + /*! partial derivative of u with respect to time at */ + inline SolScalarType u() = 0; + /*! partial derivative of u with respect to time at */ + inline SolScalarType du_dx(uint_type x_i) = 0; + /*! partial derivative of u with respect to time at */ + inline SolScalarType du_dt() = 0; + +protected: + +}; + + + +/*! provides access to the element solution vector through a memory view. */ +template +class FEVarData: public MAST::ComputeKernelBase { + +public: + + using var_scalar_type = typename MAST::DeducedScalarType::type; + using coefficient_view_type = typename ViewTraits::coefficient_view_type; + + FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm) {} + virtual ~FEVarData() {} + virtual inline void execute() override {} + inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} + inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = coeffs;} + inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} + inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) { return _du_dx(qp, comp, x_i);} + +protected: + + bool *_compute_du_dx; + const MAST::FEShapeDataBase *_fe; + const typename ViewTraits::coefficient_view_type _coeffs; + typename ViewTraits::u_view_type _u; + typename ViewTraits::du_dx_view_type _du_dx; +}; + + + + +template +class FEVarData>: +public MAST::ComputeKernelBase { + +public: + + using var_scalar_type = typename MAST::DeducedScalarType::type; + using coefficient_view_type = typename EigenFESolDataViewTraits::coefficient_view_type; + using u_view_type = typename MAST::EigenFESolDataViewTraits::u_view_type; + using du_dx_view_type = typename MAST::EigenFESolDataViewTraits::du_dx_view_type; + + + FEVarData(const std::string& nm): + MAST::ComputeKernelBase(nm), + _compute_du_dx (nullptr), + _fe (nullptr), + _coeffs (nullptr) + {} + virtual ~FEVarData() {} + + virtual inline void execute() override {} + inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} + inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = &coeffs;} + inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} + inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) + { return _du_dx(qp, _fe->spatial_dim()*comp + x_i);} + +protected: + + bool _compute_du_dx; + const MAST::FEShapeDataBase *_fe; + const coefficient_view_type *_coeffs; + u_view_type _u; + du_dx_view_type _du_dx; +}; + + + +template +class ComputeKernel: public MAST::ComputeKernelBase { + +public: + + ComputeKernel(const std::string& nm): MAST::ComputeKernelBase(nm) {} + virtual ~ComputeKernel() {} + virtual inline const ValueType& value() const = 0; +}; + + + +template +class ComputeKernelDerivative: public MAST::ComputeKernelBase { + +public: + + ComputeKernelDerivative(const std::string& nm): + MAST::ComputeKernelBase(nm), + _f (nullptr) + {} + virtual ~ComputeKernelDerivative() {} + virtual inline void set_derivative_paramter(const MAST::FunctionBase& f); + virtual inline const ValueType& value() const = 0; + +protected: + + const MAST::FunctionBase* _f; +}; + + + +template +class FEComputeKernel: public MAST::ComputeKernel { + +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; + + FEComputeKernel(const std::string& nm): MAST::ComputeKernel(nm) { } + virtual ~FEComputeKernel() { } + virtual inline const ValueType& value() const = 0; + //virtual inline const ValueType& derivative(const MAST::FunctionBase& f) const = 0; +}; + + + + +/*! Collection of compute kernels with specified dependencies and data views*/ +template +class ComputationBase: public MAST::CurrentComputation { + +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(): MAST::CurrentComputation() { } + virtual ~ComputationBase() { } + void add_compute_kernel(MAST::ComputeKernelBase& c); + /*! parses through all the compute kernels and their views to prepare a graph of dependency. */ + void prepare(); + void print_graph(); + void execute(); + +protected: + +}; + + +template +struct ElasticityComputeKernelTraits { + + using basis_scalar_type = BasisScalarType; + using nodal_scalar_type = NodalScalarType; + using sol_scalar_type = SolScalarType; + using var_scalar_type = typename MAST::DeducedScalarType; + using res_vector_type = Eigen::Matrix; + using jac_matrix_type = Eigen::Matrix; + using material_val_type = Eigen::Matrix; + using quadrature_type = MAST::libMeshQuadrature; + using fe_basis_type = MAST::libMeshFE; + using fe_derived_traits = MAST::EigenFEShapeDataViewTraits; + using fe_derived_type = MAST::FEGeometryBasisDerivedData; + using fe_var_traits = MAST::EigenFESolDataViewTraits; + using fe_var_type = typename MAST::FEVarData; +}; + + + +template +class StrainEnergy: +public FEComputeKernel { + +public: + + using value_type = typename Traits::res_vector_type; + 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; + using fe_var_traits = typename Traits::fe_var_traits; + using fe_var_type = typename Traits::fe_var_type; + using material_value_type = typename Traits::material_val_type; + + StrainEnergy(): + MAST::FEComputeKernel("strain_energy"), + _material (nullptr), + _fe_var_data (nullptr) + { } + + virtual ~StrainEnergy() { } + virtual inline void set_material(const MAST::ComputeKernel& material) + { _material = &material;} + virtual inline void set_fe_var_data(const fe_var_type& fe_data) + { _fe_var_data = &fe_data;} + virtual inline void execute() override { } + virtual inline const value_type& value() const override { } + +protected: + + + const MAST::ComputeKernel *_material; + const fe_var_type *_fe_var_data; +}; + + + + +template +class ElasticityElemOperations: public MAST::ComputationBase { + +public: + + using value_type = typename Traits::res_vector_type; + 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; + using quadrature_type = typename Traits::quadrature_type; + using fe_view_traits = typename Traits::fe_var_traits; + using fe_basis_type = typename Traits::fe_basis_type; + using fe_derived_type = typename Traits::fe_derived_type; + using fe_var_type = typename Traits::fe_var_type; + using material_value_type = typename Traits::material_val_type; + + ElasticityElemOperations(): + MAST::ComputationBase(), + _initialized (false), + _sys (nullptr), + _physics (nullptr), + _X (nullptr), + _fe_var_data (nullptr), + _strain_energy (nullptr) + { } + virtual ~ElasticityElemOperations() {} + virtual void init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& X); + virtual void clear() {} + virtual void reinit(const libMesh::Elem& e) {} + virtual void compute_residual() {} + virtual void compute_jacobian() {} + virtual void compute_complex_step_jacobian() {} + virtual void compute_auto_diff_jacobian() {} + +protected: + + bool _initialized; + const MAST::NonlinearSystem* _sys; + const MAST::PhysicsDisciplineBase* _physics; + const libMesh::NumericVector* _X; + + quadrature_type *_quadrature; + fe_basis_type *_fe_basis; + fe_derived_type *_fe_derived; + fe_var_type *_fe_var_data; + typename MAST::StrainEnergy *_strain_energy; +}; + + +template +void +MAST::ElasticityElemOperations::init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& X) { + + libmesh_assert(!_initialized); + + + // setup the volume compute kernels + _sys = &sys; + _physics = &physics; + _X = &X; + + + _quadrature = new quadrature_type("quadrature"); + _fe_basis = new fe_basis_type("fe_basis"); + _fe_basis->set_quadrature(*_quadrature); + _fe_derived = new fe_derived_type("fe_derived"); + _fe_derived->set_fe_basis(*_fe_basis); + _fe_var_data = new fe_var_type("u_vars"); + _fe_var_data->set_fe_shape_data(*_fe_derived); + _strain_energy = new MAST::StrainEnergy; + _strain_energy->set_fe_var_data(*_fe_var_data); + //_strain_energy->set_material(*_material); +} + + +} +#endif /* __mast_compute_kernel_h__ */ + From dfd38aae5b0e3a3849c35f4214394db6d30b7e61 Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Fri, 22 May 2020 22:43:10 -0500 Subject: [PATCH 2/9] -- Refactored compute kernel code to individual files. --- src/base/computation_base.h | 34 + src/base/compute_kernel.h | 781 ++---------------- src/base/compute_kernel_base.h | 53 ++ src/base/mast_data_types.h | 20 + src/base/view.h | 144 ++++ src/elasticity/elasticity_elem_operations.h | 114 +++ src/elasticity/strain_energy_compute_kernel.h | 47 ++ src/mesh/fe.h | 112 +++ src/mesh/fe_basis_derived_data.h | 102 +++ src/mesh/fe_compute_kernel.h | 25 + src/mesh/fe_geometry_basis_derived_data.h | 96 +++ src/mesh/fe_shape_data_base.h | 65 ++ src/mesh/fe_var_data.h | 80 ++ src/mesh/quadrature.h | 88 ++ 14 files changed, 1061 insertions(+), 700 deletions(-) create mode 100644 src/base/computation_base.h create mode 100644 src/base/compute_kernel_base.h create mode 100644 src/base/view.h create mode 100644 src/elasticity/elasticity_elem_operations.h create mode 100644 src/elasticity/strain_energy_compute_kernel.h create mode 100644 src/mesh/fe.h create mode 100644 src/mesh/fe_basis_derived_data.h create mode 100644 src/mesh/fe_compute_kernel.h create mode 100644 src/mesh/fe_geometry_basis_derived_data.h create mode 100644 src/mesh/fe_shape_data_base.h create mode 100644 src/mesh/fe_var_data.h create mode 100644 src/mesh/quadrature.h diff --git a/src/base/computation_base.h b/src/base/computation_base.h new file mode 100644 index 00000000..a03d29a9 --- /dev/null +++ b/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 +class ComputationBase: public MAST::CurrentComputation { + +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(): MAST::CurrentComputation() { } + virtual ~ComputationBase() { } + void add_compute_kernel(MAST::ComputeKernelBase& 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_ diff --git a/src/base/compute_kernel.h b/src/base/compute_kernel.h index 73405888..8a1c7cc6 100644 --- a/src/base/compute_kernel.h +++ b/src/base/compute_kernel.h @@ -10,491 +10,77 @@ #define _mast_compute_kernel_h_ // C++ includes -#include -#include +#include +#include // MAST includes -#include "base/mast_data_types.h" +#include "base/compute_kernel_base.h" +#include "base/view.h" -// libMesh includes -#include "libmesh/quadrature.h" -#include "libmesh/fe_base.h" -#include "libmesh/elem.h" - -typedef unsigned int uint_type; namespace MAST { +class ComputeKernelDataBase { -// Forward declerations -class FunctionBase; - - -class ComputeKernelBase { - -public: - - ComputeKernelBase(const std::string& nm): _nm(nm) {} - virtual ~ComputeKernelBase() {} - virtual inline bool depends_on(const MAST::ComputeKernelBase& d) const { return _dependency.count(&d);} - virtual inline const std::set& get_dependencies() const { return _dependency;} - virtual inline void pre_execute() {} - virtual inline void post_execute() {} - virtual inline void execute() = 0; - -protected: - - virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} - - const std::string _nm; - std::set _dependency; -}; - - -struct CurrentComputation { - - CurrentComputation(): qp(0), time (0.), dt(0.), param(nullptr) {} - - uint_type qp; - Real time; - Real dt; - /*! parameter for which sensitivity is being computed */ - FunctionBase *param; -}; - - -template -class Quadrature: public MAST::ComputeKernelBase { - -public: - - using scalar_type = ScalarType; - - Quadrature(const std::string& nm): MAST::ComputeKernelBase(nm) {} - virtual ~Quadrature() {} - virtual inline uint_type dim() const = 0; - virtual inline uint_type order() const = 0; - virtual inline uint_type n_points() const = 0; - virtual inline ScalarType qp_coord(uint_type qp, uint_type xi_i) const = 0; - virtual inline ScalarType weight(uint_type qp) const = 0; - -}; - - - -/*! This provides a */ -template -class MappedQuadrature: public MAST::Quadrature { - -public: - - using scalar_type = ScalarType; - - MappedQuadrature(const std::string& nm): MAST::Quadrature(nm), _dim(0) {} - virtual ~MappedQuadrature() {} - virtual void set_data(const uint_type dim, - const typename ViewTraits::quadrature_point_type points, - const typename ViewTraits::quadrature_point_type weights) { - _points = points; - _weights = weights; - } - virtual inline uint_type dim() const override { return _dim;} - virtual inline uint_type order() const override { /* should not be called*/ libmesh_assert(false);} - virtual inline uint_type n_points() const override { return _weights.size();} - virtual inline ScalarType qp_coord(uint_type qp, uint_type xi_i) const override { return _points(qp, xi_i);} - virtual inline ScalarType weight(uint_type qp) const { return _weights(qp);} - -protected: - - uint_type _dim; - typename ViewTraits::quadrature_point_type _points; - typename ViewTraits::quadrature_point_type _weights; -}; - - - -/*! serves as a wrapper around libMesh */ -class libMeshQuadrature: public MAST::Quadrature { - -public: - - /*! the quadrature object \p q is expected to be initialized outside of this class. */ - libMeshQuadrature(const std::string& nm): - MAST::Quadrature(nm), - _q (nullptr) - {} - virtual ~libMeshQuadrature() {} - virtual inline void execute() override { } - const libMesh::QBase& get_libmesh_object() const { return *_q;} - virtual inline uint_type dim() const override { return _q->get_dim();} - virtual inline uint_type order() const override { return _q->get_order();} - virtual inline uint_type n_points() const override { return _q->n_points();} - virtual inline scalar_type qp_coord(uint_type qp, uint_type xi_i) const override - { return _q->get_points()[qp](xi_i);} - virtual inline scalar_type weight(uint_type qp) const override { return _q->w(qp);} - -protected: - - libMesh::QBase* _q; -}; - - - -template -class FEBasis: ComputeKernelBase { - public: - using basis_scalar_type = ScalarType; - - FEBasis(const std::string& nm): - MAST::ComputeKernelBase(nm), - _quadrature (nullptr) + ComputeKernelDataBase(const std::string& nm, + const std::vector& d): + _name (nm), + _dim (d) { } - virtual ~FEBasis() {} - - virtual inline void set_quadrature(const MAST::Quadrature& q) { _quadrature = &q;} - virtual inline uint_type dim() const = 0; - virtual inline uint_type order() const = 0; - virtual inline uint_type n_basis() const = 0; - virtual inline ScalarType qp_coord(uint_type qp, uint_type x_i) const - { return _quadrature->qp_coord(qp, x_i); } - virtual inline ScalarType phi(uint_type qp, uint_type phi_i) const = 0; - virtual inline ScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; - -protected: - - const MAST::Quadrature *_quadrature; -}; - -/*! serves as a wrapper around libMesh FEBase object to provide */ -class libMeshFE: public FEBasis { - -public: - - using scalar_type = Real; - - libMeshFE(const std::string& nm): - MAST::FEBasis(nm), - _compute_dphi_dxi (false), - _fe (nullptr) + ComputeKernelDataBase(MAST::ComputeKernelDataBase& d): + _name (d._name), + _dim (d._dim) { } - virtual ~libMeshFE() { } - - inline void set_compute_dphi_dxi(bool f) { _compute_dphi_dxi = f;} - - virtual inline void reinit(const libMesh::FEBase& fe) { - _fe = &fe; - const libMesh::FEMap& m = fe.get_fe_map(); - _dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; - } - virtual inline void execute() override { } - - virtual inline uint_type dim() const override { return _fe->get_dim();} - - virtual inline uint_type order() const override {return _fe->get_order();} + virtual ~ComputeKernelDataBase() {} - virtual inline uint_type n_basis() const override { return _fe->n_shape_functions();} - - virtual inline basis_scalar_type phi(uint_type qp, uint_type phi_i) const override - {return _fe->get_phi()[phi_i][qp];} - - virtual inline basis_scalar_type - dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const override - { return (*_dphi_dxi[xi_i])[phi_i][qp];} + inline uint_type n_dim() const { return _dim.size();} protected: - - bool _compute_dphi_dxi; - const libMesh::FEBase* _fe; - std::vector>*> _dphi_dxi; -}; - - -template -struct DeducedScalarType { }; - -template <> -struct DeducedScalarType { using type = Real;}; - -template <> -struct DeducedScalarType, Real> { using type = std::complex;}; - -template <> -struct DeducedScalarType> { using type = std::complex;}; - -template <> -struct DeducedScalarType, std::complex> { using type = std::complex;}; - - -template -struct EigenFEShapeDataViewTraits { - - using xyz_view_type = Eigen::Matrix; - using detJ_view_type = Eigen::Matrix; - using detJxW_view_type = Eigen::Matrix; - using dx_dxi_view_type = Eigen::Matrix; - using dxi_dx_view_type = Eigen::Matrix; - using dphi_dx_view_type = Eigen::Matrix; - - template - void init_view(ViewType& view, uint_type n1); - - template - void init_view(ViewType& view, uint_type n1, uint_type n2); + const std::string _name; + const std::vector _dim; }; template -struct EigenFESolDataViewTraits { - - using coefficient_view_type = Eigen::Matrix; - using u_view_type = Eigen::Matrix; - using du_dx_view_type = Eigen::Matrix; - - template - void init_view(ViewType& view, uint_type n1); - - template - void init_view(ViewType& view, uint_type n1, uint_type n2); -}; - - - -template -class FEShapeDataBase: public MAST::ComputeKernelBase { - -public: - - FEShapeDataBase(const std::string& nm): - MAST::ComputeKernelBase(nm), - _compute_xyz (false), - _compute_Jac (false), - _compute_detJ (false), - _compute_JxW (false), - _compute_dphi_dx (false), - _compute_normal (false), - _spatial_dim (0), - _fe_basis (nullptr) - {} - virtual ~FEShapeDataBase() {} - inline void set_compute_xyz(bool f) { _compute_xyz = f;} - inline void set_compute_Jac(bool f) { _compute_Jac = f;} - inline void set_compute_detJ(bool f) { _compute_detJ = f;} - inline void set_compute_detJxW(bool f) { _compute_JxW = f;} - inline void set_compute_dphi_dx(bool f) { _compute_dphi_dx = f;} - inline void set_compute_normal(bool f) { _compute_normal = f;} - inline void set_fe_basis(MAST::FEBasis& basis) { _fe_basis = &basis;} - virtual inline void execute() override { } - //virtual inline void reinit_for_side(MAST::FEBasis& basis, uint_type s) = 0; - virtual inline uint_type ref_dim() const { return _fe_basis->dim();} - virtual inline uint_type spatial_dim() const { return _spatial_dim;} - virtual inline uint_type order() const { return _fe_basis->order();} - virtual inline uint_type n_basis() const { return _fe_basis->n_basis();} - virtual inline BasisScalarType phi(uint_type qp, uint_type phi_i) const - { return _fe_basis->phi(qp, phi_i);} - virtual inline BasisScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const - { return _fe_basis->dphi_dxi(qp, phi_i, xi_i);} - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const = 0; - virtual inline NodalScalarType detJ(uint_type qp) const = 0; - virtual inline NodalScalarType detJxW(uint_type qp) const = 0; - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; - -protected: - - bool _compute_xyz; - bool _compute_Jac; - bool _compute_detJ; - bool _compute_JxW; - bool _compute_dphi_dx; - bool _compute_normal; - uint_type _spatial_dim; - - MAST::FEBasis *_fe_basis; -}; - - - -/*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate - nodal locations. Typically, Lagrange shape functoins are used for this purpose */ -template -class FEGeometryBasisDerivedData: -public MAST::FEShapeDataBase { - -public: - - FEGeometryBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm) {} - virtual ~FEGeometryBasisDerivedData() {} - - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _xyz(qp, x_i);} - - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _detJ(qp);} - - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _detJxW(qp);} - - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dx_dxi(qp, x_i, xi_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dxi_dx(qp, x_i, xi_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, phi_i, x_i);} - -protected: - - typename ViewTraits::xyz_view_type _xyz; - typename ViewTraits::detJ_view_type _detJ; - typename ViewTraits::detJxW_view_type _detJxW; - typename ViewTraits::dx_dxi_view_type _dx_dxi; - typename ViewTraits::dxi_dx_view_type _dxi_dx; - typename ViewTraits::dphi_dx_view_type _dphi_dx; -}; - - - -/*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate - nodal locations. Typically, Lagrange shape functoins are used for this purpose */ -template -class FEGeometryBasisDerivedData>: -public MAST::FEShapeDataBase { +class IndexableComputeKernelData: public MAST::ComputeKernelDataBase { public: - FEGeometryBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm) {} - - virtual ~FEGeometryBasisDerivedData() {} + IndexableComputeKernelData(const std::string& nm, + const std::vector& d): + MAST::ComputeKernelDataBase(nm, d), + _view (nullptr) + { } - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _xyz(qp, x_i);} + IndexableComputeKernelData(MAST::IndexableComputeKernelData& d): + MAST::ComputeKernelDataBase(d) { } - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _detJ(qp);} - - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _detJxW(qp);} - - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} - -protected: - - typename EigenFEShapeDataViewTraits::xyz_view_type _xyz; - typename EigenFEShapeDataViewTraits::detJ_view_type _detJ; - typename EigenFEShapeDataViewTraits::detJxW_view_type _detJxW; - typename EigenFEShapeDataViewTraits::dx_dxi_view_type _dx_dxi; - typename EigenFEShapeDataViewTraits::dxi_dx_view_type _dxi_dx; - typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; -}; + virtual ~IndexableComputeKernelData() + { if (_view) delete _view; } - -/*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, - this is any non-Lagrange basis.*/ -template -class FEBasisDerivedData: public MAST::FEShapeDataBase { - -public: - - using fe_geom_derived_type = MAST::FEGeometryBasisDerivedData; - - FEBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm), - _fe_geom(nullptr) - {} + inline void set_outer_dimensions(std::vector idx, + uint_type buffer_size=0); - virtual ~FEBasisDerivedData() {} - - virtual void inline set_fe_geom_derived_data(const fe_geom_derived_type& fe_geom) - { _fe_geom = &fe_geom;} - - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _fe_geom->xyz(qp, x_i);} - - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _fe_geom->detJ(qp);} + template inline ViewType + get_local_view(const std::vector& idx) { + + // make sure this object has been initialized + libmesh_assert(_view); + + return _view->get_view_slice(idx); + } - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _fe_geom->detJxW(qp);} - - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dx_dxi(qp, x_i, xi_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dxi_dx(qp, x_i, xi_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, phi_i, x_i);} - protected: - - const fe_geom_derived_type *_fe_geom; - - typename ViewTraits::dphi_dx_view_type _dphi_dx; -}; - - - -/*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, - this is any non-Lagrange basis.*/ -template -class FEBasisDerivedData>: -public MAST::FEShapeDataBase { - -public: - - using fe_derived_data_type = MAST::FEGeometryBasisDerivedData>; - - FEBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm), - _fe_geom (nullptr) - {} - - virtual ~FEBasisDerivedData() {} - - virtual void inline set_fe_geom_derived_data(const fe_derived_data_type& fe_geom) - { _fe_geom = &fe_geom;} - - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _fe_geom->xyz(qp, x_i);} - - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _fe_geom->detJ(qp);} - - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _fe_geom->detJxW(qp);} - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dx_dxi(qp, x_i, xi_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dxi_dx(qp, x_i, xi_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} - -protected: - - const fe_derived_data_type& _fe_geom; - - typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; + MAST::View *_view; + std::vector _outer_idx_size; }; @@ -514,79 +100,12 @@ class FieldVariable { /*! partial derivative of u with respect to time at */ inline SolScalarType du_dt() = 0; -protected: - -}; - - - -/*! provides access to the element solution vector through a memory view. */ -template -class FEVarData: public MAST::ComputeKernelBase { - -public: - - using var_scalar_type = typename MAST::DeducedScalarType::type; - using coefficient_view_type = typename ViewTraits::coefficient_view_type; - - FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm) {} - virtual ~FEVarData() {} - virtual inline void execute() override {} - inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} - inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = coeffs;} - inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} - inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) { return _du_dx(qp, comp, x_i);} - -protected: - - bool *_compute_du_dx; - const MAST::FEShapeDataBase *_fe; - const typename ViewTraits::coefficient_view_type _coeffs; - typename ViewTraits::u_view_type _u; - typename ViewTraits::du_dx_view_type _du_dx; }; - -template -class FEVarData>: -public MAST::ComputeKernelBase { - -public: - - using var_scalar_type = typename MAST::DeducedScalarType::type; - using coefficient_view_type = typename EigenFESolDataViewTraits::coefficient_view_type; - using u_view_type = typename MAST::EigenFESolDataViewTraits::u_view_type; - using du_dx_view_type = typename MAST::EigenFESolDataViewTraits::du_dx_view_type; - - - FEVarData(const std::string& nm): - MAST::ComputeKernelBase(nm), - _compute_du_dx (nullptr), - _fe (nullptr), - _coeffs (nullptr) - {} - virtual ~FEVarData() {} - - virtual inline void execute() override {} - inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} - inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = &coeffs;} - inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} - inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) - { return _du_dx(qp, _fe->spatial_dim()*comp + x_i);} - -protected: - - bool _compute_du_dx; - const MAST::FEShapeDataBase *_fe; - const coefficient_view_type *_coeffs; - u_view_type _u; - du_dx_view_type _du_dx; -}; - +// Forward decleration +template class ComputeKernelDerivative; template @@ -594,211 +113,73 @@ class ComputeKernel: public MAST::ComputeKernelBase { public: - ComputeKernel(const std::string& nm): MAST::ComputeKernelBase(nm) {} - virtual ~ComputeKernel() {} - virtual inline const ValueType& value() const = 0; -}; - - - -template -class ComputeKernelDerivative: public MAST::ComputeKernelBase { - -public: + using derivative_kernel_type = MAST::ComputeKernelDerivative; - ComputeKernelDerivative(const std::string& nm): + ComputeKernel(const std::string& nm): MAST::ComputeKernelBase(nm), - _f (nullptr) + _derivative_kernel (nullptr) {} - virtual ~ComputeKernelDerivative() {} - virtual inline void set_derivative_paramter(const MAST::FunctionBase& f); - virtual inline const ValueType& value() const = 0; - -protected: - - const MAST::FunctionBase* _f; -}; - - - -template -class FEComputeKernel: public MAST::ComputeKernel { - -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; - FEComputeKernel(const std::string& nm): MAST::ComputeKernel(nm) { } - virtual ~FEComputeKernel() { } - virtual inline const ValueType& value() const = 0; - //virtual inline const ValueType& derivative(const MAST::FunctionBase& f) const = 0; -}; - - - - -/*! Collection of compute kernels with specified dependencies and data views*/ -template -class ComputationBase: public MAST::CurrentComputation { + virtual ~ComputeKernel() {} -public: + virtual inline void set_derivative_kernel(derivative_kernel_type& d) + { + libmesh_assert_msg(!_derivative_kernel, "Derivative kernel already set."); + _derivative_kernel = &d; + } - 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(): MAST::CurrentComputation() { } - virtual ~ComputationBase() { } - void add_compute_kernel(MAST::ComputeKernelBase& c); - /*! parses through all the compute kernels and their views to prepare a graph of dependency. */ - void prepare(); - void print_graph(); - void execute(); + virtual inline const derivative_kernel_type& get_derivative_kernel() const + { + libmesh_assert_msg(_derivative_kernel, "Derivative kernel not set"); + return *_derivative_kernel; + } + virtual inline ValueType value() const = 0; + protected: -}; - - -template -struct ElasticityComputeKernelTraits { - - using basis_scalar_type = BasisScalarType; - using nodal_scalar_type = NodalScalarType; - using sol_scalar_type = SolScalarType; - using var_scalar_type = typename MAST::DeducedScalarType; - using res_vector_type = Eigen::Matrix; - using jac_matrix_type = Eigen::Matrix; - using material_val_type = Eigen::Matrix; - using quadrature_type = MAST::libMeshQuadrature; - using fe_basis_type = MAST::libMeshFE; - using fe_derived_traits = MAST::EigenFEShapeDataViewTraits; - using fe_derived_type = MAST::FEGeometryBasisDerivedData; - using fe_var_traits = MAST::EigenFESolDataViewTraits; - using fe_var_type = typename MAST::FEVarData; + derivative_kernel_type *_derivative_kernel; }; -template -class StrainEnergy: -public FEComputeKernel { +template +class ComputeKernelDerivative: public MAST::ComputeKernelBase { public: - - using value_type = typename Traits::res_vector_type; - 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; - using fe_var_traits = typename Traits::fe_var_traits; - using fe_var_type = typename Traits::fe_var_type; - using material_value_type = typename Traits::material_val_type; - - StrainEnergy(): - MAST::FEComputeKernel("strain_energy"), - _material (nullptr), - _fe_var_data (nullptr) - { } - virtual ~StrainEnergy() { } - virtual inline void set_material(const MAST::ComputeKernel& material) - { _material = &material;} - virtual inline void set_fe_var_data(const fe_var_type& fe_data) - { _fe_var_data = &fe_data;} - virtual inline void execute() override { } - virtual inline const value_type& value() const override { } - -protected: + using primal_kernel_type = MAST::ComputeKernel; + + ComputeKernelDerivative (const std::string& nm): + MAST::ComputeKernelBase (nm), + _primal_kernel (nullptr), + _f (nullptr) + {} + virtual ~ComputeKernelDerivative() {} - const MAST::ComputeKernel *_material; - const fe_var_type *_fe_var_data; -}; - - - - -template -class ElasticityElemOperations: public MAST::ComputationBase { + virtual inline void set_primal_kernel(primal_kernel_type& k) + { + libmesh_assert_msg(!_primal_kernel, "Primal kernel already set."); + _primal_kernel = &k; + } -public: + virtual inline primal_kernel_type& get_primal_kernel() + { + libmesh_assert_msg(_primal_kernel, "Primal kernel not set."); + return *_primal_kernel; + } - using value_type = typename Traits::res_vector_type; - 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; - using quadrature_type = typename Traits::quadrature_type; - using fe_view_traits = typename Traits::fe_var_traits; - using fe_basis_type = typename Traits::fe_basis_type; - using fe_derived_type = typename Traits::fe_derived_type; - using fe_var_type = typename Traits::fe_var_type; - using material_value_type = typename Traits::material_val_type; - - ElasticityElemOperations(): - MAST::ComputationBase(), - _initialized (false), - _sys (nullptr), - _physics (nullptr), - _X (nullptr), - _fe_var_data (nullptr), - _strain_energy (nullptr) - { } - virtual ~ElasticityElemOperations() {} - virtual void init(const MAST::NonlinearSystem& sys, - const MAST::PhysicsDisciplineBase& physics, - const libMesh::NumericVector& X); - virtual void clear() {} - virtual void reinit(const libMesh::Elem& e) {} - virtual void compute_residual() {} - virtual void compute_jacobian() {} - virtual void compute_complex_step_jacobian() {} - virtual void compute_auto_diff_jacobian() {} + virtual inline void set_derivative_paramter(const MAST::FunctionBase& f); + virtual inline ValueType value() const = 0; protected: - bool _initialized; - const MAST::NonlinearSystem* _sys; - const MAST::PhysicsDisciplineBase* _physics; - const libMesh::NumericVector* _X; - - quadrature_type *_quadrature; - fe_basis_type *_fe_basis; - fe_derived_type *_fe_derived; - fe_var_type *_fe_var_data; - typename MAST::StrainEnergy *_strain_energy; + primal_kernel_type *_primal_kernel; + const MAST::FunctionBase *_f; }; - -template -void -MAST::ElasticityElemOperations::init(const MAST::NonlinearSystem& sys, - const MAST::PhysicsDisciplineBase& physics, - const libMesh::NumericVector& X) { - - libmesh_assert(!_initialized); - - - // setup the volume compute kernels - _sys = &sys; - _physics = &physics; - _X = &X; - - - _quadrature = new quadrature_type("quadrature"); - _fe_basis = new fe_basis_type("fe_basis"); - _fe_basis->set_quadrature(*_quadrature); - _fe_derived = new fe_derived_type("fe_derived"); - _fe_derived->set_fe_basis(*_fe_basis); - _fe_var_data = new fe_var_type("u_vars"); - _fe_var_data->set_fe_shape_data(*_fe_derived); - _strain_energy = new MAST::StrainEnergy; - _strain_energy->set_fe_var_data(*_fe_var_data); - //_strain_energy->set_material(*_material); } - -} #endif /* __mast_compute_kernel_h__ */ diff --git a/src/base/compute_kernel_base.h b/src/base/compute_kernel_base.h new file mode 100644 index 00000000..f7a6284c --- /dev/null +++ b/src/base/compute_kernel_base.h @@ -0,0 +1,53 @@ + +#ifndef _mast_compute_kernel_base_h_ +#define _mast_compute_kernel_base_h_ + +// C++ includes +#include +#include + +// MAST includes +#include "base/mast_data_types.h" + + +namespace MAST { + +// Forward declerations +class FunctionBase; + +class ComputeKernelBase { + +public: + + ComputeKernelBase(const std::string& nm): _nm(nm) {} + virtual ~ComputeKernelBase() {} + virtual inline bool depends_on(const MAST::ComputeKernelBase& d) const { return _dependency.count(&d);} + virtual inline const std::set& get_dependencies() const { return _dependency;} + virtual inline void init() {} + virtual inline void pre_execute() {} + virtual inline void post_execute() {} + virtual inline void execute() = 0; + +protected: + + virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} + + const std::string _nm; + std::set _dependency; +}; + + +struct CurrentComputation { + + CurrentComputation(): time (0.), dt(0.), param(nullptr) {} + + std::vector current_indices; + Real time; + Real dt; + /*! parameter for which sensitivity is being computed */ + FunctionBase *param; +}; + +} + +#endif // _mast_compute_kernel_base_h_ diff --git a/src/base/mast_data_types.h b/src/base/mast_data_types.h index b3afae97..4be13f05 100644 --- a/src/base/mast_data_types.h +++ b/src/base/mast_data_types.h @@ -29,6 +29,7 @@ #include "Eigen/Dense" using namespace Eigen; +typedef unsigned int uint_type; typedef libMesh::Real Real; typedef libMesh::Complex Complex; @@ -56,4 +57,23 @@ typedef libMesh::DenseMatrix DenseComplexMatrix; typedef libMesh::DenseVector DenseRealVector; typedef libMesh::DenseVector DenseComplexVector; +namespace MAST { + +template +struct DeducedScalarType { }; + +template <> +struct DeducedScalarType { using type = Real;}; + +template <> +struct DeducedScalarType, Real> { using type = std::complex;}; + +template <> +struct DeducedScalarType> { using type = std::complex;}; + +template <> +struct DeducedScalarType, std::complex> { using type = std::complex;}; + +} + #endif // __mast__data_types__ diff --git a/src/base/view.h b/src/base/view.h new file mode 100644 index 00000000..9719ca06 --- /dev/null +++ b/src/base/view.h @@ -0,0 +1,144 @@ + +#ifndef _mast_view_h_ +#define _mast_view_h_ + +// MAST includes +#include "base/mast_data_types.h" + + + +namespace MAST { + + +template +struct EigenViewType { + + using eigen_matrix_type = Eigen::Matrix; + using view_type = Eigen::Map; + using view_m1_type = typename MAST::EigenViewType::view_type; + + static inline view_type + create_view(uint_type* dim, ScalarType* array) { + uint_type + n1 = dim[0], + n2 = dim[1]; + for (unsigned int i=2; i static inline + typename MAST::EigenViewType::view_type + create_view_slice(const std::vector idx, uint_type* dim, ScalarType* array); + + template <> static inline + typename MAST::EigenViewType::view_type + create_view_slice<1>(const std::vector idx, uint_type* dim, ScalarType* array) + { + + libmesh_assert_equal_to(idx.size(), 1); + libmesh_assert_less(idx[0], dim[0]); + + uint_type + n2 = dim[1]; + + for (unsigned int i=2; i::create_view(dim+1, array+idx[0]*n2); + } +}; + + +template +struct EigenViewType { + + using eigen_matrix_type = Eigen::Matrix; + using view_type = Eigen::Map; + using view_m1_type = ScalarType; + + template static inline ScalarType + create_view_slice(const std::vector idx, uint_type* dim, ScalarType* array); + + template <> static inline ScalarType + create_view_slice<1>(const std::vector idx, uint_type* dim, ScalarType* array) { + + libmesh_assert_less(idx.size(), 1); + libmesh_assert_less(idx[0], dim[0]); + + // Dim-1 dimensional data is a scalar + return array[idx[0]]; + } + +}; + + + + +template +class View { + +public: + + View(const std::string& nm, uint_type dim, uint_type n1, ...): + _name (nm), + _dim (dim), + _dim_size (nullptr), + _array (nullptr) { + + _dim_size = new uint_type[dim]; + + va_list args; + va_start (args, n1); + + uint_type + n = 1; + + for (uint_type i=0; i inline ViewType get_view() + { return MAST::EigenViewType::create_view(_dim, _array);} + + template inline const ViewType get_view() const + { return MAST::EigenViewType::create_view(_dim, _array);} + + template inline ViewType get_view_slice(const std::vector& idx) + { return MAST::EigenViewType::create_view_slice(idx, _dim, _array);} + + template inline const ViewType get_view_slice(const std::vector& idx) const + { return MAST::EigenViewType::create_view_slice(idx, _dim, _array);} + + +protected: + + std::string _name; + uint_type _dim; + uint_type *_dim_size; + ScalarType *_array; +}; + + + +} + +#endif // _mast_view_h_ diff --git a/src/elasticity/elasticity_elem_operations.h b/src/elasticity/elasticity_elem_operations.h new file mode 100644 index 00000000..f2e74d36 --- /dev/null +++ b/src/elasticity/elasticity_elem_operations.h @@ -0,0 +1,114 @@ + +#ifndef _mast_elasticity_elem_operations_h_ +#define _mast_elasticity_elem_operations_h_ + + +// MAST includes +#include "base/computation_base.h" +#include "mesh/quadrature.h" +#include "mesh/fe.h" +#include "mesh/fe_geometry_basis_derived_data.h" +#include "mesh/fe_var_data.h" +#include "elasticity/strain_energy_compute_kernel.h" + +namespace MAST { + +template +struct ElasticityComputeKernelTraits { + + using basis_scalar_type = BasisScalarType; + using nodal_scalar_type = NodalScalarType; + using sol_scalar_type = SolScalarType; + using var_scalar_type = typename MAST::DeducedScalarType; + using res_vector_type = Eigen::Matrix; + using jac_matrix_type = Eigen::Matrix; + using material_val_type = Eigen::Matrix; + using quadrature_type = MAST::libMeshQuadrature; + using fe_basis_type = MAST::libMeshFE; + using fe_derived_traits = MAST::EigenFEShapeDataViewTraits; + using fe_derived_type = MAST::FEGeometryBasisDerivedData; + using fe_var_traits = MAST::EigenFESolDataViewTraits; + using fe_var_type = typename MAST::FEVarData; +}; + + +template +class ElasticityElemOperations: public MAST::ComputationBase { + +public: + + using value_type = typename Traits::res_vector_type; + 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; + using quadrature_type = typename Traits::quadrature_type; + using fe_view_traits = typename Traits::fe_var_traits; + using fe_basis_type = typename Traits::fe_basis_type; + using fe_derived_type = typename Traits::fe_derived_type; + using fe_var_type = typename Traits::fe_var_type; + using material_value_type = typename Traits::material_val_type; + + ElasticityElemOperations(): + MAST::ComputationBase(), + _initialized (false), + _sys (nullptr), + _physics (nullptr), + _X (nullptr), + _fe_var_data (nullptr), + _strain_energy (nullptr) + { } + virtual ~ElasticityElemOperations() {} + virtual void init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& X); + virtual void clear() {} + virtual void reinit(const libMesh::Elem& e) {} + virtual void compute_residual() {} + virtual void compute_jacobian() {} + virtual void compute_complex_step_jacobian() {} + virtual void compute_auto_diff_jacobian() {} + +protected: + + bool _initialized; + const MAST::NonlinearSystem* _sys; + const MAST::PhysicsDisciplineBase* _physics; + const libMesh::NumericVector* _X; + + quadrature_type *_quadrature; + fe_basis_type *_fe_basis; + fe_derived_type *_fe_derived; + fe_var_type *_fe_var_data; + typename MAST::StrainEnergy *_strain_energy; +}; + + +template +void +MAST::ElasticityElemOperations::init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& X) { + + libmesh_assert(!_initialized); + + + // setup the volume compute kernels + _sys = &sys; + _physics = &physics; + _X = &X; + + + _quadrature = new quadrature_type("quadrature"); + _fe_basis = new fe_basis_type("fe_basis"); + _fe_basis->set_quadrature(*_quadrature); + _fe_derived = new fe_derived_type("fe_derived"); + _fe_derived->set_fe_basis(*_fe_basis); + _fe_var_data = new fe_var_type("u_vars"); + _fe_var_data->set_fe_shape_data(*_fe_derived); + _strain_energy = new MAST::StrainEnergy; + _strain_energy->set_fe_var_data(*_fe_var_data); + //_strain_energy->set_material(*_material); +} +} + +#endif // _mast_elasticity_elem_operations_h_ diff --git a/src/elasticity/strain_energy_compute_kernel.h b/src/elasticity/strain_energy_compute_kernel.h new file mode 100644 index 00000000..2992b2f0 --- /dev/null +++ b/src/elasticity/strain_energy_compute_kernel.h @@ -0,0 +1,47 @@ + +#ifndef _mast_strain_energy_compute_kernel_h_ +#define _mast_strain_energy_compute_kernel_h_ + +// MAST includes +#include "mesh/fe_compute_kernel.h" + +namespace MAST { + +template +class StrainEnergy: +public FEComputeKernel { + +public: + + using value_type = typename Traits::res_vector_type; + 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; + using fe_var_traits = typename Traits::fe_var_traits; + using fe_var_type = typename Traits::fe_var_type; + using material_value_type = typename Traits::material_val_type; + + StrainEnergy(): + MAST::FEComputeKernel("strain_energy"), + _material (nullptr), + _fe_var_data (nullptr) + { } + + virtual ~StrainEnergy() { } + virtual inline void set_material(const MAST::ComputeKernel& material) + { _material = &material;} + virtual inline void set_fe_var_data(const fe_var_type& fe_data) + { _fe_var_data = &fe_data;} + virtual inline void execute() override { } + virtual inline value_type value() const override { } + +protected: + + + const MAST::ComputeKernel *_material; + const fe_var_type *_fe_var_data; +}; + +} + +#endif // _mast_strain_energy_compute_kernel_h_ diff --git a/src/mesh/fe.h b/src/mesh/fe.h new file mode 100644 index 00000000..3b644f1a --- /dev/null +++ b/src/mesh/fe.h @@ -0,0 +1,112 @@ + +#ifndef _mast_febasis_h_ +#define _mast_febasis_h_ + +// MAST includes +#include "base/compute_kernel.h" +#include "mesh/quadrature.h" + +namespace MAST { + +template +struct EigenFEShapeDataViewTraits { + + using xyz_view_type = Eigen::Matrix; + using detJ_view_type = Eigen::Matrix; + using detJxW_view_type = Eigen::Matrix; + using dx_dxi_view_type = Eigen::Matrix; + using dxi_dx_view_type = Eigen::Matrix; + using dphi_dx_view_type = Eigen::Matrix; +}; + + + +template +struct EigenFESolDataViewTraits { + + using coefficient_view_type = Eigen::Matrix; + using u_view_type = Eigen::Matrix; + using du_dx_view_type = Eigen::Matrix; +}; + + + + +template +class FEBasis: ComputeKernelBase { + +public: + + using basis_scalar_type = ScalarType; + + FEBasis(const std::string& nm): + MAST::ComputeKernelBase(nm), + _quadrature (nullptr) + { } + virtual ~FEBasis() {} + + virtual inline void set_quadrature(const MAST::Quadrature& q) { _quadrature = &q;} + virtual inline uint_type dim() const = 0; + virtual inline uint_type order() const = 0; + virtual inline uint_type n_basis() const = 0; + virtual inline ScalarType qp_coord(uint_type qp, uint_type x_i) const + { return _quadrature->qp_coord(qp, x_i); } + virtual inline ScalarType phi(uint_type qp, uint_type phi_i) const = 0; + virtual inline ScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + +protected: + + const MAST::Quadrature *_quadrature; +}; + + +/*! serves as a wrapper around libMesh FEBase object to provide */ +class libMeshFE: public FEBasis { + +public: + + using scalar_type = Real; + + libMeshFE(const std::string& nm): + MAST::FEBasis(nm), + _compute_dphi_dxi (false), + _fe (nullptr) + { } + virtual ~libMeshFE() { } + + inline void set_compute_dphi_dxi(bool f) { _compute_dphi_dxi = f;} + + virtual inline void reinit(const libMesh::FEBase& fe) { + _fe = &fe; + const libMesh::FEMap& m = fe.get_fe_map(); + _dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; + } + + virtual inline void execute() override { } + + virtual inline uint_type dim() const override { return _fe->get_dim();} + + virtual inline uint_type order() const override {return _fe->get_order();} + + virtual inline uint_type n_basis() const override { return _fe->n_shape_functions();} + + virtual inline basis_scalar_type phi(uint_type qp, uint_type phi_i) const override + {return _fe->get_phi()[phi_i][qp];} + + virtual inline basis_scalar_type + dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const override + { return (*_dphi_dxi[xi_i])[phi_i][qp];} + +protected: + + bool _compute_dphi_dxi; + const libMesh::FEBase* _fe; + std::vector>*> _dphi_dxi; +}; + + + + +} + +#endif // _mast_febasis_h_ diff --git a/src/mesh/fe_basis_derived_data.h b/src/mesh/fe_basis_derived_data.h new file mode 100644 index 00000000..c4c2884d --- /dev/null +++ b/src/mesh/fe_basis_derived_data.h @@ -0,0 +1,102 @@ + +#ifndef _mast_fe_basis_derived_data_h_ +#define _mast_fe_basis_derived_data_h_ + +// MAST includes +#include "mesh/fe_shape_data_base.h" + +namespace MAST { + +/*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, + this is any non-Lagrange basis.*/ +template +class FEBasisDerivedData: public MAST::FEShapeDataBase { + +public: + + using fe_geom_derived_type = MAST::FEGeometryBasisDerivedData; + + FEBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm), + _fe_geom(nullptr) + {} + + virtual ~FEBasisDerivedData() {} + + virtual void inline set_fe_geom_derived_data(const fe_geom_derived_type& fe_geom) + { _fe_geom = &fe_geom;} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _fe_geom->xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _fe_geom->detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _fe_geom->detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dx_dxi(qp, x_i, xi_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dxi_dx(qp, x_i, xi_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, phi_i, x_i);} + +protected: + + const fe_geom_derived_type *_fe_geom; + + typename ViewTraits::dphi_dx_view_type _dphi_dx; +}; + + + +/*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, + this is any non-Lagrange basis.*/ +template +class FEBasisDerivedData>: +public MAST::FEShapeDataBase { + +public: + + using fe_derived_data_type = MAST::FEGeometryBasisDerivedData>; + + FEBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm), + _fe_geom (nullptr) + {} + + virtual ~FEBasisDerivedData() {} + + virtual void inline set_fe_geom_derived_data(const fe_derived_data_type& fe_geom) + { _fe_geom = &fe_geom;} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _fe_geom->xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _fe_geom->detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _fe_geom->detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dx_dxi(qp, x_i, xi_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _fe_geom->dxi_dx(qp, x_i, xi_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} + +protected: + + const fe_derived_data_type& _fe_geom; + + typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; +}; +} + +#endif // _mast_fe_basis_derived_data_h_ diff --git a/src/mesh/fe_compute_kernel.h b/src/mesh/fe_compute_kernel.h new file mode 100644 index 00000000..fbef430e --- /dev/null +++ b/src/mesh/fe_compute_kernel.h @@ -0,0 +1,25 @@ + +#ifndef _mast_fe_compute_kernel_h_ +#define _mast_fe_compute_kernel_h_ + +// MAST includes +#include "base/compute_kernel.h" + +namespace MAST { + +template +class FEComputeKernel: public MAST::ComputeKernel< ValueType> { + +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; + + FEComputeKernel(const std::string& nm): MAST::ComputeKernel(nm) { } + virtual ~FEComputeKernel() { } +}; + +} + +#endif // _mast_fe_compute_kernel_h_ diff --git a/src/mesh/fe_geometry_basis_derived_data.h b/src/mesh/fe_geometry_basis_derived_data.h new file mode 100644 index 00000000..fd1b2b33 --- /dev/null +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -0,0 +1,96 @@ + +#ifndef _mast_fe_geometry_basis_derived_data_h_ +#define _mast_fe_geometry_basis_derived_data_h_ + +// MAST includes +#include "mesh/fe_shape_data_base.h" + + +namespace MAST { + +/*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate + nodal locations. Typically, Lagrange shape functoins are used for this purpose */ +template +class FEGeometryBasisDerivedData: +public MAST::FEShapeDataBase { + +public: + + FEGeometryBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm) {} + virtual ~FEGeometryBasisDerivedData() {} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dx_dxi(qp, x_i, xi_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dxi_dx(qp, x_i, xi_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, phi_i, x_i);} + +protected: + + typename ViewTraits::xyz_view_type _xyz; + typename ViewTraits::detJ_view_type _detJ; + typename ViewTraits::detJxW_view_type _detJxW; + typename ViewTraits::dx_dxi_view_type _dx_dxi; + typename ViewTraits::dxi_dx_view_type _dxi_dx; + typename ViewTraits::dphi_dx_view_type _dphi_dx; +}; + + + +/*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate + nodal locations. Typically, Lagrange shape functoins are used for this purpose */ +template +class FEGeometryBasisDerivedData>: +public MAST::FEShapeDataBase { + +public: + + FEGeometryBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm) {} + + virtual ~FEGeometryBasisDerivedData() {} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override + { return _xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const override + { return _detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { return _detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override + { return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i);} + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override + { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} + +protected: + + typename EigenFEShapeDataViewTraits::xyz_view_type _xyz; + typename EigenFEShapeDataViewTraits::detJ_view_type _detJ; + typename EigenFEShapeDataViewTraits::detJxW_view_type _detJxW; + typename EigenFEShapeDataViewTraits::dx_dxi_view_type _dx_dxi; + typename EigenFEShapeDataViewTraits::dxi_dx_view_type _dxi_dx; + typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; +}; + +} + +#endif // _mast_fe_geometry_basis_derived_data_h_ diff --git a/src/mesh/fe_shape_data_base.h b/src/mesh/fe_shape_data_base.h new file mode 100644 index 00000000..ee7245b4 --- /dev/null +++ b/src/mesh/fe_shape_data_base.h @@ -0,0 +1,65 @@ + +#ifndef _mast_fe_shape_data_base_h_ +#define _mast_fe_shape_data_base_h_ + +// MAST includes +#include "mesh/fe.h" + +namespace MAST { + +template +class FEShapeDataBase: public MAST::ComputeKernelBase { + +public: + + FEShapeDataBase(const std::string& nm): + MAST::ComputeKernelBase(nm), + _compute_xyz (false), + _compute_Jac (false), + _compute_detJ (false), + _compute_JxW (false), + _compute_dphi_dx (false), + _compute_normal (false), + _spatial_dim (0), + _fe_basis (nullptr) + {} + virtual ~FEShapeDataBase() {} + inline void set_compute_xyz(bool f) { _compute_xyz = f;} + inline void set_compute_Jac(bool f) { _compute_Jac = f;} + inline void set_compute_detJ(bool f) { _compute_detJ = f;} + inline void set_compute_detJxW(bool f) { _compute_JxW = f;} + inline void set_compute_dphi_dx(bool f) { _compute_dphi_dx = f;} + inline void set_compute_normal(bool f) { _compute_normal = f;} + inline void set_fe_basis(MAST::FEBasis& basis) { _fe_basis = &basis;} + virtual inline void execute() override { } + //virtual inline void reinit_for_side(MAST::FEBasis& basis, uint_type s) = 0; + virtual inline uint_type ref_dim() const { return _fe_basis->dim();} + virtual inline uint_type spatial_dim() const { return _spatial_dim;} + virtual inline uint_type order() const { return _fe_basis->order();} + virtual inline uint_type n_basis() const { return _fe_basis->n_basis();} + virtual inline BasisScalarType phi(uint_type qp, uint_type phi_i) const + { return _fe_basis->phi(qp, phi_i);} + virtual inline BasisScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const + { return _fe_basis->dphi_dxi(qp, phi_i, xi_i);} + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const = 0; + virtual inline NodalScalarType detJ(uint_type qp) const = 0; + virtual inline NodalScalarType detJxW(uint_type qp) const = 0; + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + +protected: + + bool _compute_xyz; + bool _compute_Jac; + bool _compute_detJ; + bool _compute_JxW; + bool _compute_dphi_dx; + bool _compute_normal; + uint_type _spatial_dim; + + MAST::FEBasis *_fe_basis; +}; + +} +#endif // _mast_fe_shape_data_base_h_ diff --git a/src/mesh/fe_var_data.h b/src/mesh/fe_var_data.h new file mode 100644 index 00000000..e7ac7999 --- /dev/null +++ b/src/mesh/fe_var_data.h @@ -0,0 +1,80 @@ + +#ifndef _mast_fe_var_data_h_ +#define _mast_fe_var_data_h_ + +// MAST includes +#include "base/compute_kernel_base.h" + + +namespace MAST { + +/*! provides access to the element solution vector through a memory view. */ +template +class FEVarData: public MAST::ComputeKernelBase { + +public: + + using var_scalar_type = typename MAST::DeducedScalarType::type; + using coefficient_view_type = typename ViewTraits::coefficient_view_type; + + FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm) {} + virtual ~FEVarData() {} + virtual inline void execute() override {} + inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} + inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = coeffs;} + inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} + inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) { return _du_dx(qp, comp, x_i);} + +protected: + + bool *_compute_du_dx; + const MAST::FEShapeDataBase *_fe; + const typename ViewTraits::coefficient_view_type _coeffs; + typename ViewTraits::u_view_type _u; + typename ViewTraits::du_dx_view_type _du_dx; +}; + + + + +template +class FEVarData>: +public MAST::ComputeKernelBase { + +public: + + using var_scalar_type = typename MAST::DeducedScalarType::type; + using coefficient_view_type = typename EigenFESolDataViewTraits::coefficient_view_type; + using u_view_type = typename MAST::EigenFESolDataViewTraits::u_view_type; + using du_dx_view_type = typename MAST::EigenFESolDataViewTraits::du_dx_view_type; + + + FEVarData(const std::string& nm): + MAST::ComputeKernelBase(nm), + _compute_du_dx (nullptr), + _fe (nullptr), + _coeffs (nullptr) + {} + virtual ~FEVarData() {} + + virtual inline void execute() override {} + inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} + inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = &coeffs;} + inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} + inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) + { return _du_dx(qp, _fe->spatial_dim()*comp + x_i);} + +protected: + + bool _compute_du_dx; + const MAST::FEShapeDataBase *_fe; + const coefficient_view_type *_coeffs; + u_view_type _u; + du_dx_view_type _du_dx; +}; + +} + +#endif // _mast_fe_var_data_h_ diff --git a/src/mesh/quadrature.h b/src/mesh/quadrature.h new file mode 100644 index 00000000..a286a246 --- /dev/null +++ b/src/mesh/quadrature.h @@ -0,0 +1,88 @@ + + +#ifndef _mast_quadrature_h_ +#define _mast_quadrature_h_ + +// MAST includes +#include "base/compute_kernel.h" + + +namespace MAST { + +template +class Quadrature: public MAST::ComputeKernelBase { + +public: + + using scalar_type = ScalarType; + + Quadrature(const std::string& nm): MAST::ComputeKernelBase(nm) {} + virtual ~Quadrature() {} + virtual inline uint_type dim() const = 0; + virtual inline uint_type order() const = 0; + virtual inline uint_type n_points() const = 0; + virtual inline ScalarType qp_coord(uint_type qp, uint_type xi_i) const = 0; + virtual inline ScalarType weight(uint_type qp) const = 0; + +}; + + + +/*! This provides a */ +template +class MappedQuadrature: public MAST::Quadrature { + +public: + + using scalar_type = ScalarType; + + MappedQuadrature(const std::string& nm): MAST::Quadrature(nm), _dim(0) {} + virtual ~MappedQuadrature() {} + virtual void set_data(const uint_type dim, + const typename ViewTraits::quadrature_point_type points, + const typename ViewTraits::quadrature_point_type weights) { + _points = points; + _weights = weights; + } + virtual inline uint_type dim() const override { return _dim;} + virtual inline uint_type order() const override { /* should not be called*/ libmesh_assert(false);} + virtual inline uint_type n_points() const override { return _weights.size();} + virtual inline ScalarType qp_coord(uint_type qp, uint_type xi_i) const override { return _points(qp, xi_i);} + virtual inline ScalarType weight(uint_type qp) const { return _weights(qp);} + +protected: + + uint_type _dim; + typename ViewTraits::quadrature_point_type _points; + typename ViewTraits::quadrature_point_type _weights; +}; + + + +/*! serves as a wrapper around libMesh */ +class libMeshQuadrature: public MAST::Quadrature { + +public: + + /*! the quadrature object \p q is expected to be initialized outside of this class. */ + libMeshQuadrature(const std::string& nm): + MAST::Quadrature(nm), + _q (nullptr) + {} + virtual ~libMeshQuadrature() {} + virtual inline void execute() override { } + const libMesh::QBase& get_libmesh_object() const { return *_q;} + virtual inline uint_type dim() const override { return _q->get_dim();} + virtual inline uint_type order() const override { return _q->get_order();} + virtual inline uint_type n_points() const override { return _q->n_points();} + virtual inline scalar_type qp_coord(uint_type qp, uint_type xi_i) const override + { return _q->get_points()[qp](xi_i);} + virtual inline scalar_type weight(uint_type qp) const override { return _q->w(qp);} + +protected: + + libMesh::QBase* _q; +}; +} + +#endif //_mast_quadrature_h_ From a0e432c821d57b491f454e02101d5f5797dc7806 Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Fri, 22 May 2020 22:43:41 -0500 Subject: [PATCH 3/9] Material property compute kernels --- src/materials/hookean_material.h | 31 ++++++++ src/materials/isotropic_hookean_material.h | 85 ++++++++++++++++++++++ 2 files changed, 116 insertions(+) create mode 100644 src/materials/hookean_material.h create mode 100644 src/materials/isotropic_hookean_material.h diff --git a/src/materials/hookean_material.h b/src/materials/hookean_material.h new file mode 100644 index 00000000..154d62ff --- /dev/null +++ b/src/materials/hookean_material.h @@ -0,0 +1,31 @@ + + +#ifndef _mast_hookean_material_h_ +#define _mast_hookean_material_h_ + +// MAST includes +#include "mesh/fe_compute_kernel.h" + +namespace MAST { + +template +class HookeanMaterial: +public ComputeKernel { + +public: + + using value_type = typename Traits::material_value_type; + + HookeanMaterial(const std::string& nm): + MAST::ComputeKernel(nm) + { } + + virtual ~HookeanMaterial() { } + +protected: + +}; + +} + +#endif // _mast_hookean_material_h_ diff --git a/src/materials/isotropic_hookean_material.h b/src/materials/isotropic_hookean_material.h new file mode 100644 index 00000000..675ad040 --- /dev/null +++ b/src/materials/isotropic_hookean_material.h @@ -0,0 +1,85 @@ + + +#ifndef _mast_isotropic_hookean_material_h_ +#define _mast_isotropic_hookean_material_h_ + +// MAST includes +#include "mesh/fe_compute_kernel.h" + +namespace MAST { + +template +class IsotropicHookeanMaterial: +public ComputeKernel { + +public: + + using scalar_type = typename Traits::var_scalar_type; + using value_type = typename Traits::material_value_type; + + IsotropicHookeanMaterial(const std::string& nm): + MAST::ComputeKernel(nm), + _E (nullptr), + _nu (nullptr) + { } + + virtual ~IsotropicHookeanMaterial() { } + + virtual inline void set_properties(const MAST::ComputeKernel& E, + const MAST::ComputeKernel& nu) { + + libmesh_assert_msg(!_E && !_nu, "Property kernels should be cleared before setting."); + + _E = &E; + _nu = ν + } + + virtual inline void execute() override { } + + + virtual inline void execute_derivative(const std::vector& idx) override { + + libmesh_assert_msg(_E && _nu, "Property kernels should be provided before execute()."); + + scalar_type + E = _E->value(), + nu = _nu->value(), + G = E/2./(1.+nu); + + MAST::init_view(_matrix, 6, 6, true); + + m(0, 0) = m(1, 1) = m(2, 2) = E*(1.-nu)/(1.-nu-2.*nu*nu); + m(0, 1) = m(0, 2) = m(1, 0) = m(1, 2) = m(2, 0) = m(2, 1) = E*nu/(1.-nu-2.*nu*nu); + m(3, 3) = m(4, 4) = m(5, 5) = E/2./(1.+nu); + } + + virtual inline value_type value() const override { + + libmesh_assert_msg(_E && _nu, "Property kernels should be provided before execute()."); + + scalar_type + E = _E->value(), + nu = _nu->value(), + G = E/2./(1.+nu); + + value_type matrix = + + MAST::init_view(_matrix, 6, 6, true); + + m(0, 0) = m(1, 1) = m(2, 2) = E*(1.-nu)/(1.-nu-2.*nu*nu); + m(0, 1) = m(0, 2) = m(1, 0) = m(1, 2) = m(2, 0) = m(2, 1) = E*nu/(1.-nu-2.*nu*nu); + m(3, 3) = m(4, 4) = m(5, 5) = E/2./(1.+nu); + } + + virtual inline value_type value_derivative() const override { return _matrix;} + +protected: + + const MAST::ComputeKernel *_E; + const MAST::ComputeKernel *_nu; +}; + +} + +#endif // _mast_isotropic_hookean_material_h_ + From d2a5608099df8c73476a587a4d3361f0cb22b5cf Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Fri, 22 May 2020 23:39:23 -0500 Subject: [PATCH 4/9] Added context_type to compute_kernels --- src/base/computation_base.h | 8 ++-- src/base/compute_kernel.h | 22 ++++----- src/base/compute_kernel_base.h | 18 ++++---- src/elasticity/elasticity_elem_operations.h | 33 ++++++++----- src/elasticity/strain_energy_compute_kernel.h | 15 +++--- src/mesh/fe.h | 24 ++++++---- src/mesh/fe_compute_kernel.h | 7 +-- src/mesh/fe_geometry_basis_derived_data.h | 14 +++--- src/mesh/fe_shape_data_base.h | 13 +++--- src/mesh/fe_var_data.h | 46 ++++++++++--------- src/mesh/quadrature.h | 23 ++++++---- 11 files changed, 124 insertions(+), 99 deletions(-) diff --git a/src/base/computation_base.h b/src/base/computation_base.h index a03d29a9..7bf2f096 100644 --- a/src/base/computation_base.h +++ b/src/base/computation_base.h @@ -8,8 +8,8 @@ namespace MAST { /*! Collection of compute kernels with specified dependencies and data views*/ -template -class ComputationBase: public MAST::CurrentComputation { +template +class ComputationBase { public: @@ -17,9 +17,9 @@ class ComputationBase: public MAST::CurrentComputation { using nodal_scalar_type = typename Traits::nodal_scalar_type; using sol_scalar_type = typename Traits::sol_scalar_type; - ComputationBase(): MAST::CurrentComputation() { } + ComputationBase() { } virtual ~ComputationBase() { } - void add_compute_kernel(MAST::ComputeKernelBase& c); + void add_compute_kernel(MAST::ComputeKernelBase& c); /*! parses through all the compute kernels and their views to prepare a graph of dependency. */ void prepare(); void print_graph(); diff --git a/src/base/compute_kernel.h b/src/base/compute_kernel.h index 8a1c7cc6..af7737c4 100644 --- a/src/base/compute_kernel.h +++ b/src/base/compute_kernel.h @@ -105,18 +105,18 @@ class FieldVariable { // Forward decleration -template class ComputeKernelDerivative; +template class ComputeKernelDerivative; -template -class ComputeKernel: public MAST::ComputeKernelBase { +template +class ComputeKernel: public MAST::ComputeKernelBase { public: - using derivative_kernel_type = MAST::ComputeKernelDerivative; + using derivative_kernel_type = MAST::ComputeKernelDerivative; ComputeKernel(const std::string& nm): - MAST::ComputeKernelBase(nm), + MAST::ComputeKernelBase(nm), _derivative_kernel (nullptr) {} @@ -134,7 +134,7 @@ class ComputeKernel: public MAST::ComputeKernelBase { return *_derivative_kernel; } - virtual inline ValueType value() const = 0; + virtual inline ValueType value(ContextType& c) const = 0; protected: @@ -143,15 +143,15 @@ class ComputeKernel: public MAST::ComputeKernelBase { -template -class ComputeKernelDerivative: public MAST::ComputeKernelBase { +template +class ComputeKernelDerivative: public MAST::ComputeKernelBase { public: - using primal_kernel_type = MAST::ComputeKernel; + using primal_kernel_type = MAST::ComputeKernel; ComputeKernelDerivative (const std::string& nm): - MAST::ComputeKernelBase (nm), + MAST::ComputeKernelBase (nm), _primal_kernel (nullptr), _f (nullptr) {} @@ -171,7 +171,7 @@ class ComputeKernelDerivative: public MAST::ComputeKernelBase { } virtual inline void set_derivative_paramter(const MAST::FunctionBase& f); - virtual inline ValueType value() const = 0; + virtual inline ValueType value(ContextType& c) const = 0; protected: diff --git a/src/base/compute_kernel_base.h b/src/base/compute_kernel_base.h index f7a6284c..74b9ce2d 100644 --- a/src/base/compute_kernel_base.h +++ b/src/base/compute_kernel_base.h @@ -15,25 +15,27 @@ namespace MAST { // Forward declerations class FunctionBase; +template class ComputeKernelBase { public: ComputeKernelBase(const std::string& nm): _nm(nm) {} virtual ~ComputeKernelBase() {} - virtual inline bool depends_on(const MAST::ComputeKernelBase& d) const { return _dependency.count(&d);} - virtual inline const std::set& get_dependencies() const { return _dependency;} - virtual inline void init() {} - virtual inline void pre_execute() {} - virtual inline void post_execute() {} - virtual inline void execute() = 0; + virtual inline bool depends_on(const MAST::ComputeKernelBase& d) const + { return _dependency.count(&d);} + virtual inline const std::set*>& get_dependencies() const + { return _dependency;} + virtual inline void pre_execute(ContextType& c) {} + virtual inline void post_execute(ContextType& c) {} + virtual inline void execute(ContextType& c) = 0; protected: - virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} + virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} const std::string _nm; - std::set _dependency; + std::set*> _dependency; }; diff --git a/src/elasticity/elasticity_elem_operations.h b/src/elasticity/elasticity_elem_operations.h index f2e74d36..d89dbef3 100644 --- a/src/elasticity/elasticity_elem_operations.h +++ b/src/elasticity/elasticity_elem_operations.h @@ -13,9 +13,10 @@ namespace MAST { -template +template struct ElasticityComputeKernelTraits { + using context_type = ContextType; using basis_scalar_type = BasisScalarType; using nodal_scalar_type = NodalScalarType; using sol_scalar_type = SolScalarType; @@ -23,20 +24,27 @@ struct ElasticityComputeKernelTraits { using res_vector_type = Eigen::Matrix; using jac_matrix_type = Eigen::Matrix; using material_val_type = Eigen::Matrix; - using quadrature_type = MAST::libMeshQuadrature; - using fe_basis_type = MAST::libMeshFE; + using quadrature_type = MAST::libMeshQuadrature; + using fe_basis_type = MAST::libMeshFE; using fe_derived_traits = MAST::EigenFEShapeDataViewTraits; - using fe_derived_type = MAST::FEGeometryBasisDerivedData; + using fe_derived_type = MAST::FEGeometryBasisDerivedData; using fe_var_traits = MAST::EigenFESolDataViewTraits; - using fe_var_type = typename MAST::FEVarData; + using fe_var_type = typename MAST::FEVarData; +}; + + +class ElasticityElemOperationsContext { + }; template -class ElasticityElemOperations: public MAST::ComputationBase { +class ElasticityElemOperations: +public MAST::ComputationBase { public: + using context_type = typename Traits::context_type; using value_type = typename Traits::res_vector_type; using basis_scalar_type = typename Traits::basis_scalar_type; using nodal_scalar_type = typename Traits::nodal_scalar_type; @@ -49,7 +57,7 @@ class ElasticityElemOperations: public MAST::ComputationBase { using material_value_type = typename Traits::material_val_type; ElasticityElemOperations(): - MAST::ComputationBase(), + MAST::ComputationBase(), _initialized (false), _sys (nullptr), _physics (nullptr), @@ -79,15 +87,16 @@ class ElasticityElemOperations: public MAST::ComputationBase { fe_basis_type *_fe_basis; fe_derived_type *_fe_derived; fe_var_type *_fe_var_data; - typename MAST::StrainEnergy *_strain_energy; + typename MAST::StrainEnergy *_strain_energy; }; template void -MAST::ElasticityElemOperations::init(const MAST::NonlinearSystem& sys, - const MAST::PhysicsDisciplineBase& physics, - const libMesh::NumericVector& X) { +MAST::ElasticityElemOperations:: +init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& X) { libmesh_assert(!_initialized); @@ -105,7 +114,7 @@ MAST::ElasticityElemOperations::init(const MAST::NonlinearSystem& sys, _fe_derived->set_fe_basis(*_fe_basis); _fe_var_data = new fe_var_type("u_vars"); _fe_var_data->set_fe_shape_data(*_fe_derived); - _strain_energy = new MAST::StrainEnergy; + _strain_energy = new MAST::StrainEnergy; _strain_energy->set_fe_var_data(*_fe_var_data); //_strain_energy->set_material(*_material); } diff --git a/src/elasticity/strain_energy_compute_kernel.h b/src/elasticity/strain_energy_compute_kernel.h index 2992b2f0..647b8d3f 100644 --- a/src/elasticity/strain_energy_compute_kernel.h +++ b/src/elasticity/strain_energy_compute_kernel.h @@ -7,9 +7,9 @@ namespace MAST { -template +template class StrainEnergy: -public FEComputeKernel { +public FEComputeKernel { public: @@ -22,23 +22,24 @@ public FEComputeKernel { using material_value_type = typename Traits::material_val_type; StrainEnergy(): - MAST::FEComputeKernel("strain_energy"), + MAST::FEComputeKernel("strain_energy"), _material (nullptr), _fe_var_data (nullptr) { } virtual ~StrainEnergy() { } - virtual inline void set_material(const MAST::ComputeKernel& material) + virtual inline void + set_material(const MAST::ComputeKernel& material) { _material = &material;} virtual inline void set_fe_var_data(const fe_var_type& fe_data) { _fe_var_data = &fe_data;} - virtual inline void execute() override { } - virtual inline value_type value() const override { } + virtual inline void execute(ContextType& c) override { } + virtual inline value_type value(ContextType& c) const override { } protected: - const MAST::ComputeKernel *_material; + const MAST::ComputeKernel *_material; const fe_var_type *_fe_var_data; }; diff --git a/src/mesh/fe.h b/src/mesh/fe.h index 3b644f1a..10c2716b 100644 --- a/src/mesh/fe.h +++ b/src/mesh/fe.h @@ -32,20 +32,21 @@ struct EigenFESolDataViewTraits { -template -class FEBasis: ComputeKernelBase { +template +class FEBasis: ComputeKernelBase { public: using basis_scalar_type = ScalarType; FEBasis(const std::string& nm): - MAST::ComputeKernelBase(nm), + MAST::ComputeKernelBase(nm), _quadrature (nullptr) { } virtual ~FEBasis() {} - virtual inline void set_quadrature(const MAST::Quadrature& q) { _quadrature = &q;} + virtual inline void set_quadrature(const MAST::Quadrature& q) + { _quadrature = &q;} virtual inline uint_type dim() const = 0; virtual inline uint_type order() const = 0; virtual inline uint_type n_basis() const = 0; @@ -56,19 +57,21 @@ class FEBasis: ComputeKernelBase { protected: - const MAST::Quadrature *_quadrature; + const MAST::Quadrature *_quadrature; }; /*! serves as a wrapper around libMesh FEBase object to provide */ -class libMeshFE: public FEBasis { +template +class libMeshFE: public FEBasis { public: - using scalar_type = Real; + using scalar_type = Real; + using basis_scalar_type = typename MAST::FEBasis::basis_scalar_type; libMeshFE(const std::string& nm): - MAST::FEBasis(nm), + MAST::FEBasis(nm), _compute_dphi_dxi (false), _fe (nullptr) { } @@ -79,10 +82,11 @@ class libMeshFE: public FEBasis { virtual inline void reinit(const libMesh::FEBase& fe) { _fe = &fe; const libMesh::FEMap& m = fe.get_fe_map(); - _dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; + this->_dphi_dxi = + {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; } - virtual inline void execute() override { } + virtual inline void execute(ContextType& c) override { } virtual inline uint_type dim() const override { return _fe->get_dim();} diff --git a/src/mesh/fe_compute_kernel.h b/src/mesh/fe_compute_kernel.h index fbef430e..95becb0f 100644 --- a/src/mesh/fe_compute_kernel.h +++ b/src/mesh/fe_compute_kernel.h @@ -7,8 +7,8 @@ namespace MAST { -template -class FEComputeKernel: public MAST::ComputeKernel< ValueType> { +template +class FEComputeKernel: public MAST::ComputeKernel { public: @@ -16,7 +16,8 @@ class FEComputeKernel: public MAST::ComputeKernel< ValueType> { using nodal_scalar_type = typename Traits::nodal_scalar_type; using sol_scalar_type = typename Traits::sol_scalar_type; - FEComputeKernel(const std::string& nm): MAST::ComputeKernel(nm) { } + FEComputeKernel(const std::string& nm): + MAST::ComputeKernel(nm) { } virtual ~FEComputeKernel() { } }; diff --git a/src/mesh/fe_geometry_basis_derived_data.h b/src/mesh/fe_geometry_basis_derived_data.h index fd1b2b33..67c49705 100644 --- a/src/mesh/fe_geometry_basis_derived_data.h +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -10,14 +10,14 @@ namespace MAST { /*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate nodal locations. Typically, Lagrange shape functoins are used for this purpose */ -template +template class FEGeometryBasisDerivedData: -public MAST::FEShapeDataBase { +public MAST::FEShapeDataBase { public: FEGeometryBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm) {} + MAST::FEShapeDataBase(nm) {} virtual ~FEGeometryBasisDerivedData() {} virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override @@ -52,14 +52,14 @@ public MAST::FEShapeDataBase { /*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate nodal locations. Typically, Lagrange shape functoins are used for this purpose */ -template -class FEGeometryBasisDerivedData>: -public MAST::FEShapeDataBase { +template +class FEGeometryBasisDerivedData, ContextType>: +public MAST::FEShapeDataBase { public: FEGeometryBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm) {} + MAST::FEShapeDataBase(nm) {} virtual ~FEGeometryBasisDerivedData() {} diff --git a/src/mesh/fe_shape_data_base.h b/src/mesh/fe_shape_data_base.h index ee7245b4..244ca124 100644 --- a/src/mesh/fe_shape_data_base.h +++ b/src/mesh/fe_shape_data_base.h @@ -7,13 +7,13 @@ namespace MAST { -template -class FEShapeDataBase: public MAST::ComputeKernelBase { +template +class FEShapeDataBase: public MAST::ComputeKernelBase { public: FEShapeDataBase(const std::string& nm): - MAST::ComputeKernelBase(nm), + MAST::ComputeKernelBase(nm), _compute_xyz (false), _compute_Jac (false), _compute_detJ (false), @@ -30,8 +30,9 @@ class FEShapeDataBase: public MAST::ComputeKernelBase { inline void set_compute_detJxW(bool f) { _compute_JxW = f;} inline void set_compute_dphi_dx(bool f) { _compute_dphi_dx = f;} inline void set_compute_normal(bool f) { _compute_normal = f;} - inline void set_fe_basis(MAST::FEBasis& basis) { _fe_basis = &basis;} - virtual inline void execute() override { } + inline void set_fe_basis(MAST::FEBasis& basis) + { _fe_basis = &basis;} + virtual inline void execute(ContextType& c) override { } //virtual inline void reinit_for_side(MAST::FEBasis& basis, uint_type s) = 0; virtual inline uint_type ref_dim() const { return _fe_basis->dim();} virtual inline uint_type spatial_dim() const { return _spatial_dim;} @@ -58,7 +59,7 @@ class FEShapeDataBase: public MAST::ComputeKernelBase { bool _compute_normal; uint_type _spatial_dim; - MAST::FEBasis *_fe_basis; + MAST::FEBasis *_fe_basis; }; } diff --git a/src/mesh/fe_var_data.h b/src/mesh/fe_var_data.h index e7ac7999..db37dbc6 100644 --- a/src/mesh/fe_var_data.h +++ b/src/mesh/fe_var_data.h @@ -9,38 +9,40 @@ namespace MAST { /*! provides access to the element solution vector through a memory view. */ -template -class FEVarData: public MAST::ComputeKernelBase { +template +class FEVarData: public MAST::ComputeKernelBase { public: using var_scalar_type = typename MAST::DeducedScalarType::type; using coefficient_view_type = typename ViewTraits::coefficient_view_type; + using fe_shape_data_type = typename MAST::FEShapeDataBase; - FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm) {} + FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm) {} virtual ~FEVarData() {} - virtual inline void execute() override {} + virtual inline void execute(ContextType& c) override {} inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} + inline void + set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = coeffs;} inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) { return _du_dx(qp, comp, x_i);} protected: - bool *_compute_du_dx; - const MAST::FEShapeDataBase *_fe; - const typename ViewTraits::coefficient_view_type _coeffs; - typename ViewTraits::u_view_type _u; - typename ViewTraits::du_dx_view_type _du_dx; + bool *_compute_du_dx; + const fe_shape_data_type *_fe; + const typename ViewTraits::coefficient_view_type _coeffs; + typename ViewTraits::u_view_type _u; + typename ViewTraits::du_dx_view_type _du_dx; }; -template -class FEVarData>: -public MAST::ComputeKernelBase { +template +class FEVarData, ContextType>: +public MAST::ComputeKernelBase { public: @@ -48,19 +50,19 @@ public MAST::ComputeKernelBase { using coefficient_view_type = typename EigenFESolDataViewTraits::coefficient_view_type; using u_view_type = typename MAST::EigenFESolDataViewTraits::u_view_type; using du_dx_view_type = typename MAST::EigenFESolDataViewTraits::du_dx_view_type; - + using fe_shape_data_type = typename MAST::FEShapeDataBase; FEVarData(const std::string& nm): - MAST::ComputeKernelBase(nm), + MAST::ComputeKernelBase(nm), _compute_du_dx (nullptr), _fe (nullptr), _coeffs (nullptr) {} virtual ~FEVarData() {} - virtual inline void execute() override {} + virtual inline void execute(ContextType& c) override {} inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - inline void set_fe_shape_data(const MAST::FEShapeDataBase& fe) { _fe = &fe;} + inline void set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = &coeffs;} inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) @@ -68,11 +70,11 @@ public MAST::ComputeKernelBase { protected: - bool _compute_du_dx; - const MAST::FEShapeDataBase *_fe; - const coefficient_view_type *_coeffs; - u_view_type _u; - du_dx_view_type _du_dx; + bool _compute_du_dx; + const fe_shape_data_type *_fe; + const coefficient_view_type *_coeffs; + u_view_type _u; + du_dx_view_type _du_dx; }; } diff --git a/src/mesh/quadrature.h b/src/mesh/quadrature.h index a286a246..6b837f73 100644 --- a/src/mesh/quadrature.h +++ b/src/mesh/quadrature.h @@ -9,14 +9,15 @@ namespace MAST { -template -class Quadrature: public MAST::ComputeKernelBase { +template +class Quadrature: public MAST::ComputeKernelBase { public: using scalar_type = ScalarType; - Quadrature(const std::string& nm): MAST::ComputeKernelBase(nm) {} + Quadrature(const std::string& nm): + MAST::ComputeKernelBase(nm) {} virtual ~Quadrature() {} virtual inline uint_type dim() const = 0; virtual inline uint_type order() const = 0; @@ -29,14 +30,15 @@ class Quadrature: public MAST::ComputeKernelBase { /*! This provides a */ -template -class MappedQuadrature: public MAST::Quadrature { +template +class MappedQuadrature: public MAST::Quadrature { public: using scalar_type = ScalarType; - MappedQuadrature(const std::string& nm): MAST::Quadrature(nm), _dim(0) {} + MappedQuadrature(const std::string& nm): + MAST::Quadrature(nm), _dim(0) {} virtual ~MappedQuadrature() {} virtual void set_data(const uint_type dim, const typename ViewTraits::quadrature_point_type points, @@ -60,17 +62,20 @@ class MappedQuadrature: public MAST::Quadrature { /*! serves as a wrapper around libMesh */ -class libMeshQuadrature: public MAST::Quadrature { +template +class libMeshQuadrature: public MAST::Quadrature { public: + using scalar_type = typename MAST::Quadrature::scalar_type; + /*! the quadrature object \p q is expected to be initialized outside of this class. */ libMeshQuadrature(const std::string& nm): - MAST::Quadrature(nm), + MAST::Quadrature(nm), _q (nullptr) {} virtual ~libMeshQuadrature() {} - virtual inline void execute() override { } + virtual inline void execute(ContextType& c) override { } const libMesh::QBase& get_libmesh_object() const { return *_q;} virtual inline uint_type dim() const override { return _q->get_dim();} virtual inline uint_type order() const override { return _q->get_order();} From 84cafe8ca14fde349fc9702d57ecd19c33252142 Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Fri, 22 May 2020 23:41:26 -0500 Subject: [PATCH 5/9] -- compute kernel instantiations added to example 6 to debug compilation. --- examples/structural/example_6/example_6.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp index 14b1c3d0..60d121c0 100644 --- a/examples/structural/example_6/example_6.cpp +++ b/examples/structural/example_6/example_6.cpp @@ -51,7 +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 "base/compute_kernel.h" +#include "elasticity/elasticity_elem_operations.h" +#include "base/view.h" // libMesh includes @@ -1449,10 +1450,10 @@ int main(int argc, char* argv[]) { optimizer->optimize(); } - MAST::ElasticityElemOperations> elem_ops; - MAST::ElasticityElemOperations> elem_ops_complex_shape; - MAST::ElasticityElemOperations> elem_ops_complex_sol; - MAST::ElasticityElemOperations> elem_ops_complex_sol2; + MAST::ElasticityElemOperations> elem_ops; + MAST::ElasticityElemOperations> elem_ops_complex_shape; + MAST::ElasticityElemOperations> elem_ops_complex_sol; + MAST::ElasticityElemOperations> elem_ops_complex_sol2; // END_TRANSLATE return 0; From b26e10f040961e48539a937978dc99ba705d8fb9 Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Wed, 10 Jun 2020 00:45:51 -0500 Subject: [PATCH 6/9] -- Cleaned up the template code. -- Multiple bug-fixes. -- Templatized FEOperatorMatrix -- Added classes for strain energy, isotropic material property. --- examples/structural/example_6/example_6.cpp | 8 +- src/base/CMakeLists.txt | 1 + src/base/compute_kernel.h | 124 +++--- src/base/compute_kernel_base.h | 12 +- src/base/mast_data_types.h | 37 ++ src/base/view.h | 156 ++++++-- src/elasticity/elasticity_elem_operations.h | 71 ++-- src/elasticity/strain_energy_compute_kernel.h | 219 ++++++++-- src/materials/isotropic_hookean_material.h | 36 +- src/mesh/fe.h | 103 +++-- src/mesh/fe_basis_derived_data.h | 100 +++-- src/mesh/fe_compute_kernel.h | 6 +- src/mesh/fe_geometry_basis_derived_data.h | 78 ++-- src/mesh/fe_shape_data_base.h | 120 +++++- src/mesh/fe_var_data.h | 60 ++- src/mesh/quadrature.h | 67 +++- src/numerics/fem_operator_matrix.cpp | 13 - src/numerics/fem_operator_matrix.h | 378 ++++++++++-------- 18 files changed, 993 insertions(+), 596 deletions(-) diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp index 60d121c0..9e017ae3 100644 --- a/examples/structural/example_6/example_6.cpp +++ b/examples/structural/example_6/example_6.cpp @@ -1450,10 +1450,10 @@ int main(int argc, char* argv[]) { optimizer->optimize(); } - MAST::ElasticityElemOperations> elem_ops; - MAST::ElasticityElemOperations> elem_ops_complex_shape; - MAST::ElasticityElemOperations> elem_ops_complex_sol; - MAST::ElasticityElemOperations> elem_ops_complex_sol2; + MAST::ElasticityElemOperations> elem_ops; + MAST::ElasticityElemOperations> elem_ops_complex_shape; + MAST::ElasticityElemOperations> elem_ops_complex_sol; + MAST::ElasticityElemOperations> elem_ops_complex_sol2; // END_TRANSLATE return 0; diff --git a/src/base/CMakeLists.txt b/src/base/CMakeLists.txt index 75ad7ea4..6d039eb5 100644 --- a/src/base/CMakeLists.txt +++ b/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 diff --git a/src/base/compute_kernel.h b/src/base/compute_kernel.h index af7737c4..f548f430 100644 --- a/src/base/compute_kernel.h +++ b/src/base/compute_kernel.h @@ -20,103 +20,80 @@ namespace MAST { -class ComputeKernelDataBase { -public: - - ComputeKernelDataBase(const std::string& nm, - const std::vector& d): - _name (nm), - _dim (d) - { } +template +struct IsScalarType { + using type = std::false_type; +}; - ComputeKernelDataBase(MAST::ComputeKernelDataBase& d): - _name (d._name), - _dim (d._dim) - { } - virtual ~ComputeKernelDataBase() {} - - inline uint_type n_dim() const { return _dim.size();} +template <> +struct IsScalarType { + using type = std::true_type; +}; -protected: - const std::string _name; - const std::vector _dim; +template <> +struct IsScalarType> { + using type = std::true_type; }; -template -class IndexableComputeKernelData: public MAST::ComputeKernelDataBase { - -public: - - IndexableComputeKernelData(const std::string& nm, - const std::vector& d): - MAST::ComputeKernelDataBase(nm, d), - _view (nullptr) - { } - - IndexableComputeKernelData(MAST::IndexableComputeKernelData& d): - MAST::ComputeKernelDataBase(d) { } - - virtual ~IndexableComputeKernelData() - { if (_view) delete _view; } +template ::type> +struct KernelReturnType { +}; - inline void set_outer_dimensions(std::vector idx, - uint_type buffer_size=0); - - template inline ViewType - get_local_view(const std::vector& idx) { - - // make sure this object has been initialized - libmesh_assert(_view); - - return _view->get_view_slice(idx); - } + +template +struct KernelReturnType { -protected: + using type = KernelValueType; +}; + - MAST::View *_view; - std::vector _outer_idx_size; +template +struct KernelReturnType { + + using type = KernelValueType; }; -/*! provides access to the element solution vector through a memory view. */ -template -class FieldVariable { - -public: +template +struct KernelReturnType { + + using type = Eigen::Map; +}; - FieldVariable(); - virtual ~FieldVariable(); - /*! partial derivative of u with respect to time at */ - inline SolScalarType u() = 0; - /*! partial derivative of u with respect to time at */ - inline SolScalarType du_dx(uint_type x_i) = 0; - /*! partial derivative of u with respect to time at */ - inline SolScalarType du_dt() = 0; +template +struct KernelReturnType { + + using type = KernelValueType&; }; +template +KernelViewType build_kernel_view(KernelType& k, ContextType& c) { } + + + // Forward decleration -template class ComputeKernelDerivative; +template class ComputeKernelDerivative; -template +template class ComputeKernel: public MAST::ComputeKernelBase { public: - using derivative_kernel_type = MAST::ComputeKernelDerivative; + using derivative_kernel_type = MAST::ComputeKernelDerivative; - ComputeKernel(const std::string& nm): - MAST::ComputeKernelBase(nm), + ComputeKernel(const std::string& nm, const bool executable): + MAST::ComputeKernelBase(nm, executable), _derivative_kernel (nullptr) {} @@ -133,24 +110,22 @@ class ComputeKernel: public MAST::ComputeKernelBase { libmesh_assert_msg(_derivative_kernel, "Derivative kernel not set"); return *_derivative_kernel; } - - virtual inline ValueType value(ContextType& c) const = 0; - + protected: - derivative_kernel_type *_derivative_kernel; + derivative_kernel_type *_derivative_kernel; }; -template +template class ComputeKernelDerivative: public MAST::ComputeKernelBase { public: - using primal_kernel_type = MAST::ComputeKernel; + using primal_kernel_type = MAST::ComputeKernel; - ComputeKernelDerivative (const std::string& nm): + ComputeKernelDerivative (const std::string& nm, const bool executable): MAST::ComputeKernelBase (nm), _primal_kernel (nullptr), _f (nullptr) @@ -171,7 +146,6 @@ class ComputeKernelDerivative: public MAST::ComputeKernelBase { } virtual inline void set_derivative_paramter(const MAST::FunctionBase& f); - virtual inline ValueType value(ContextType& c) const = 0; protected: @@ -179,6 +153,8 @@ class ComputeKernelDerivative: public MAST::ComputeKernelBase { const MAST::FunctionBase *_f; }; + + } #endif /* __mast_compute_kernel_h__ */ diff --git a/src/base/compute_kernel_base.h b/src/base/compute_kernel_base.h index 74b9ce2d..6e2bca20 100644 --- a/src/base/compute_kernel_base.h +++ b/src/base/compute_kernel_base.h @@ -20,21 +20,25 @@ class ComputeKernelBase { public: - ComputeKernelBase(const std::string& nm): _nm(nm) {} + 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& d) const { return _dependency.count(&d);} virtual inline const std::set*>& get_dependencies() const { return _dependency;} - virtual inline void pre_execute(ContextType& c) {} - virtual inline void post_execute(ContextType& c) {} - virtual inline void execute(ContextType& c) = 0; protected: virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} const std::string _nm; + const bool _executable; std::set*> _dependency; }; diff --git a/src/base/mast_data_types.h b/src/base/mast_data_types.h index 4be13f05..ae020c41 100644 --- a/src/base/mast_data_types.h +++ b/src/base/mast_data_types.h @@ -30,6 +30,7 @@ using namespace Eigen; typedef unsigned int uint_type; +typedef int int_type; typedef libMesh::Real Real; typedef libMesh::Complex Complex; @@ -57,6 +58,42 @@ typedef libMesh::DenseMatrix DenseComplexMatrix; typedef libMesh::DenseVector DenseRealVector; typedef libMesh::DenseVector DenseComplexVector; +template +struct EigenVector +{ + using type = Eigen::Matrix; + using map_type = Eigen::Map; +}; + +template +struct EigenRowVector +{ + using type = Eigen::Matrix; + using map_type = Eigen::Map; +}; + + +template +struct EigenMatrix +{ + using type = Eigen::Matrix; + using map_type = Eigen::Map; +}; + + +struct EigenTraits { + + template + using matrix_type = typename EigenMatrix::type; + template + using matrix_map_type = typename EigenMatrix::map_type; + template + using vector_type = typename EigenVector::type; + template + using vector_map_type = typename EigenVector::map_type; +}; + + namespace MAST { template diff --git a/src/base/view.h b/src/base/view.h index 9719ca06..e2b5e051 100644 --- a/src/base/view.h +++ b/src/base/view.h @@ -2,6 +2,9 @@ #ifndef _mast_view_h_ #define _mast_view_h_ +// C++ includes +#include + // MAST includes #include "base/mast_data_types.h" @@ -76,45 +79,99 @@ struct EigenViewType { -template -class View { +class ViewBase { + public: - - View(const std::string& nm, uint_type dim, uint_type n1, ...): + ViewBase(const std::string& nm, const std::vector& sizes): _name (nm), - _dim (dim), - _dim_size (nullptr), - _array (nullptr) { + _dim (sizes.size()), + _size (sizes) { - _dim_size = new uint_type[dim]; + // this should be a vector or higher dimensional data + libmesh_assert_greater(_dim, 0); + + _slice_size.resize(_dim+1, 1); + + for (uint_type i=_dim-1; i>=0; i++) + _slice_size[i] = _slice_size[i+1]*_size[i]; + } + + virtual ~ViewBase() {} - va_list args; - va_start (args, n1); + inline const std::string& name() const { return _name;} - uint_type - n = 1; + virtual inline void set_outer_sizes(const std::vector& d) + { + libmesh_assert_equal_to(_outside_size.size(), 0); + _outside_size = d; + _total_size.resize(d.size()+_size.size(), 0); + for (uint_type i=0; i& indices) { + + libmesh_assert_equal_to(indices.size(), _dim); - for (uint_type i=0; i _size; + std::vector _slice_size; + std::vector _outside_size; + std::vector _total_size; +}; + + +template +class View: public MAST::ViewBase { + +public: + + View(const std::string& nm, const std::vector& sizes): + MAST::ViewBase(nm, sizes), + _array (nullptr) { - virtual ~View() - { - delete _dim_size; - delete _array; + } + + virtual ~View() { delete _array;} + + virtual inline void init () override { + + libmesh_assert_msg(!_array, "View array already initialized."); + + _array = new ScalarType[_slice_size[0]]; + + for (uint_type i=0; i<_slice_size[0]; i++) + _array[i] = ScalarType{}; } + inline ScalarType& operator() (const std::vector& indices) { + + libmesh_assert_equal_to(indices.size(), _dim); + return _array[array_index(indices)]; + } + + inline const ScalarType operator() (const std::vector& indices) const { + + libmesh_assert_equal_to(indices.size(), _dim); + return _array[array_index(indices)]; + } template inline ViewType get_view() { return MAST::EigenViewType::create_view(_dim, _array);} @@ -128,16 +185,55 @@ class View { template inline const ViewType get_view_slice(const std::vector& idx) const { return MAST::EigenViewType::create_view_slice(idx, _dim, _array);} - protected: - std::string _name; - uint_type _dim; - uint_type *_dim_size; - ScalarType *_array; + ScalarType *_array; }; +class ComputationData { + +public: + + ComputationData(){} + + virtual ~ComputationData() {} + + inline void add_data(MAST::ViewBase& v) { + + const std::string& nm = v.name(); + + libmesh_assert_msg(_data_map.find(nm) == _data_map.end(), "Duplicate data name: "+nm); + + _data_map[nm] = &v; + } + + template MAST::View& get_data(const std::string& nm) { + + std::map::iterator + it = _data_map.find(nm), + end = _data_map.end(); + + libmesh_assert_msg( it != end, "Invalid data name: "+nm); + return static_cast&>(*it->second); + } + + virtual inline void set_outer_sizes_and_initialize(const std::vector& d) + { + std::map::iterator + it = _data_map.begin(), + end = _data_map.end(); + + for ( ; it != end; it++) { + it->second->set_outer_sizes(d); + it->second->init(); + } + } + +protected: + + std::map _data_map; +}; } diff --git a/src/elasticity/elasticity_elem_operations.h b/src/elasticity/elasticity_elem_operations.h index d89dbef3..2d6b7a43 100644 --- a/src/elasticity/elasticity_elem_operations.h +++ b/src/elasticity/elasticity_elem_operations.h @@ -13,28 +13,43 @@ namespace MAST { +template +class ConstantSectionMaterial { +public: + + ConstantSectionMaterial(){} + virtual ~ConstantSectionMaterial(){} + virtual inline void value(const ContextType& c, typename EigenMatrix::type& m) const { + libmesh_assert_equal_to(m.rows(), 3); + libmesh_assert_equal_to(m.cols(), 3); + m(0,0) = 1.; m(0, 1) = .5; + m(1,0) = .5; m(1, 1) = 1.; + m(2,2) = 0.75; + } +}; + + template -struct ElasticityComputeKernelTraits { +struct ElasticityTraits { - using context_type = ContextType; - using basis_scalar_type = BasisScalarType; - using nodal_scalar_type = NodalScalarType; - using sol_scalar_type = SolScalarType; - using var_scalar_type = typename MAST::DeducedScalarType; - using res_vector_type = Eigen::Matrix; - using jac_matrix_type = Eigen::Matrix; - using material_val_type = Eigen::Matrix; - using quadrature_type = MAST::libMeshQuadrature; - using fe_basis_type = MAST::libMeshFE; - using fe_derived_traits = MAST::EigenFEShapeDataViewTraits; - using fe_derived_type = MAST::FEGeometryBasisDerivedData; - using fe_var_traits = MAST::EigenFESolDataViewTraits; - using fe_var_type = typename MAST::FEVarData; + using context_type = ContextType; + using basis_scalar_type = BasisScalarType; + using nodal_scalar_type = NodalScalarType; + using sol_scalar_type = SolScalarType; + using var_scalar_type = typename MAST::DeducedScalarType::type; + using view_traits = EigenTraits; + using section_property_type = MAST::ConstantSectionMaterial; + using quadrature_type = MAST::libMeshQuadrature; + using fe_basis_type = MAST::libMeshFE; + using fe_derived_type = MAST::FEGeometryBasisDerivedData; + using fe_var_type = MAST::FEVarData; }; -class ElasticityElemOperationsContext { +struct ElasticityElemOperationsContext { + ElasticityElemOperationsContext(): qp(0) {} + uint_type qp; }; @@ -44,17 +59,15 @@ public MAST::ComputationBase { public: - using context_type = typename Traits::context_type; - using value_type = typename Traits::res_vector_type; - 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; - using quadrature_type = typename Traits::quadrature_type; - using fe_view_traits = typename Traits::fe_var_traits; - using fe_basis_type = typename Traits::fe_basis_type; - using fe_derived_type = typename Traits::fe_derived_type; - using fe_var_type = typename Traits::fe_var_type; - using material_value_type = typename Traits::material_val_type; + using context_type = typename Traits::context_type; + 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; + using quadrature_type = typename Traits::quadrature_type; + using fe_basis_type = typename Traits::fe_basis_type; + using fe_derived_type = typename Traits::fe_derived_type; + using fe_var_type = typename Traits::fe_var_type; + using section_property_type = typename Traits::section_property_type; ElasticityElemOperations(): MAST::ComputationBase(), @@ -87,7 +100,7 @@ public MAST::ComputationBase { fe_basis_type *_fe_basis; fe_derived_type *_fe_derived; fe_var_type *_fe_var_data; - typename MAST::StrainEnergy *_strain_energy; + typename MAST::Element2D::LinearContinuumStrainEnergy *_strain_energy; }; @@ -114,7 +127,7 @@ init(const MAST::NonlinearSystem& sys, _fe_derived->set_fe_basis(*_fe_basis); _fe_var_data = new fe_var_type("u_vars"); _fe_var_data->set_fe_shape_data(*_fe_derived); - _strain_energy = new MAST::StrainEnergy; + _strain_energy = new MAST::Element2D::LinearContinuumStrainEnergy; _strain_energy->set_fe_var_data(*_fe_var_data); //_strain_energy->set_material(*_material); } diff --git a/src/elasticity/strain_energy_compute_kernel.h b/src/elasticity/strain_energy_compute_kernel.h index 647b8d3f..888416fb 100644 --- a/src/elasticity/strain_energy_compute_kernel.h +++ b/src/elasticity/strain_energy_compute_kernel.h @@ -4,45 +4,216 @@ // MAST includes #include "mesh/fe_compute_kernel.h" +#include "numerics/fem_operator_matrix.h" + namespace MAST { -template -class StrainEnergy: -public FEComputeKernel { +namespace Element2D { + +template +void linear_continuum_strain(const FEVarType& fe_var, + const uint_type qp, + typename EigenVector::type& epsilon, + MAST::FEMOperatorMatrixBase& Bmat_lin) { + + epsilon.setZero(); + + const typename FEVarType::fe_shape_data_type + &fe = fe_var.get_fe_shape_data(); + + // make sure all matrices are the right size + libmesh_assert_equal_to(epsilon.size(), 3); + libmesh_assert_equal_to(Bmat_lin.m(), 3); + libmesh_assert_equal_to(Bmat_lin.n(), 2*fe.n_basis()); + + + // linear strain operator + Bmat_lin.set_shape_function(0, 0, fe.dphi_dx(qp, 0)); // epsilon_xx = du/dx + Bmat_lin.set_shape_function(2, 1, fe.dphi_dx(qp, 0)); // gamma_xy = dv/dx + ... + + // linear strain operator + Bmat_lin.set_shape_function(1, 1, fe.dphi_dx(qp, 1)); // epsilon_yy = dv/dy + Bmat_lin.set_shape_function(2, 0, fe.dphi_dx(qp, 1)); // gamma_xy = du/dy + ... + + epsilon(0) = fe_var.du_dx(qp, 0, 0); // du/dx + epsilon(1) = fe_var.du_dx(qp, 1, 1); // dv/dy + epsilon(2) = fe_var.du_dx(qp, 0, 1) + fe_var.du_dx(qp, 1, 0); // du/dy + dv/dx +} + + +template +void green_lagrange_strain_operator(const FEVarType& fe_var, + const ContextType& qp, + bool if_nonlinear, + EigenMatrix& E, + EigenMatrix& F, + EigenVector& epsilon, + EigenMatrix& mat_x, + EigenMatrix& mat_y, + MAST::FEMOperatorMatrixBase& Bmat_lin, + MAST::FEMOperatorMatrixBase& Bmat_nl_x, + MAST::FEMOperatorMatrixBase& Bmat_nl_y, + MAST::FEMOperatorMatrixBase& Bmat_nl_u, + MAST::FEMOperatorMatrixBase& Bmat_nl_v) { + + + epsilon.setZero(); + mat_x.setZero(); + mat_y.setZero(); + + const typename FEVarType::fe_shape_data_type + &fe = fe_var.get_fe_shape_object(); + + // make sure all matrices are the right size + libmesh_assert_equal_to(epsilon.size(), 3); + libmesh_assert_equal_to(mat_x.rows(), 3); + libmesh_assert_equal_to(mat_x.cols(), 2); + libmesh_assert_equal_to(mat_y.rows(), 3); + libmesh_assert_equal_to(mat_y.cols(), 2); + libmesh_assert_equal_to(Bmat_lin.m(), 3); + libmesh_assert_equal_to(Bmat_lin.n(), 2*fe.n_basis()); + libmesh_assert_equal_to(Bmat_nl_x.m(), 2); + libmesh_assert_equal_to(Bmat_nl_x.n(), 2*fe.n_basis()); + libmesh_assert_equal_to(Bmat_nl_y.m(), 2); + libmesh_assert_equal_to(Bmat_nl_y.n(), 2*fe.n_basis()); + libmesh_assert_equal_to(F.cols(), 2); + libmesh_assert_equal_to(F.rows(), 2); + libmesh_assert_equal_to(E.cols(), 2); + libmesh_assert_equal_to(E.rows(), 2); + + + // now set the shape function values + Bmat_lin.set_shape_function(0, 0, fe.dphi_dx(qp, 0)); // epsilon_xx = du/dx + Bmat_lin.set_shape_function(2, 1, fe.dphi_dx(qp, 0)); // gamma_xy = dv/dx + ... + + // nonlinear strain operator in x + Bmat_nl_x.set_shape_function(0, 0, fe.dphi_dx(qp, 0)); // du/dx + Bmat_nl_x.set_shape_function(1, 1, fe.dphi_dx(qp, 0)); // dv/dx + + // nonlinear strain operator in u + Bmat_nl_u.set_shape_function(0, 0, fe.dphi_dx(qp, 0)); // du/dx + Bmat_nl_v.set_shape_function(0, 1, fe.dphi_dx(qp, 0)); // dv/dx + + // dN/dy + Bmat_lin.set_shape_function(1, 1, fe.dphi_dx(qp, 1)); // epsilon_yy = dv/dy + Bmat_lin.set_shape_function(2, 0, fe.dphi_dx(qp, 1)); // gamma_xy = du/dy + ... + + // nonlinear strain operator in y + Bmat_nl_y.set_shape_function(0, 0, fe.dphi_dx(qp, 1)); // du/dy + Bmat_nl_y.set_shape_function(1, 1, fe.dphi_dx(qp, 1)); // dv/dy + + // nonlinear strain operator in v + Bmat_nl_u.set_shape_function(1, 0, fe.dphi_dx(qp, 1)); // du/dy + Bmat_nl_v.set_shape_function(1, 1, fe.dphi_dx(qp, 1)); // dv/dy + + // prepare the deformation gradient matrix + F.row(0) = fe_var.du_dx(qp, 0); + F.row(1) = fe_var.du_dx(qp, 1); + + // this calculates the Green-Lagrange strain in the reference config + E = 0.5*(F + F.transpose() + F.transpose() * F); + + // now, add this to the strain vector + epsilon(0) = E(0,0); + epsilon(1) = E(1,1); + epsilon(2) = E(0,1) + E(1,0); + + // now initialize the matrices with strain components + // that multiply the Bmat_nl terms + mat_x(0, 0) = fe_var.du_dx(qp, 0, 0); + mat_x(0, 1) = fe_var.du_dx(qp, 1, 0); + mat_x(2, 0) = fe_var.du_dx(qp, 0, 1); + mat_x(2, 1) = fe_var.du_dx(qp, 1, 1); + + mat_y(1, 0) = fe_var.du_dx(qp, 0, 1); + mat_y(1, 1) = fe_var.du_dx(qp, 1, 1); + mat_y(2, 0) = fe_var.du_dx(qp, 0, 0); + mat_y(2, 1) = fe_var.du_dx(qp, 1, 0); +} -public: - using value_type = typename Traits::res_vector_type; - 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; - using fe_var_traits = typename Traits::fe_var_traits; - using fe_var_type = typename Traits::fe_var_type; - using material_value_type = typename Traits::material_val_type; - - StrainEnergy(): - MAST::FEComputeKernel("strain_energy"), - _material (nullptr), +template +class LinearContinuumStrainEnergy: +public ComputeKernel { + +public: + + using nodal_scalar_type = typename Traits::nodal_scalar_type; + using var_scalar_type = typename Traits::var_scalar_type; + using vector_type = typename Traits::view_traits::template vector_type; + using matrix_type = typename Traits::view_traits::template matrix_type; + using fe_var_type = typename Traits::fe_var_type; + using section_property_type = typename Traits::section_property_type; + + LinearContinuumStrainEnergy(): + MAST::ComputeKernel("strain_energy", false), + _property (nullptr), _fe_var_data (nullptr) { } - virtual ~StrainEnergy() { } + virtual ~LinearContinuumStrainEnergy() { } + virtual inline void - set_material(const MAST::ComputeKernel& material) - { _material = &material;} + set_section_property(const section_property_type& p) { + + libmesh_assert_msg(!_property, "Property already initialized."); + + _property = &p; + } + virtual inline void set_fe_var_data(const fe_var_type& fe_data) - { _fe_var_data = &fe_data;} - virtual inline void execute(ContextType& c) override { } - virtual inline value_type value(ContextType& c) const override { } + { + libmesh_assert_msg(!_fe_var_data, "FE data already initialized."); + _fe_var_data = &fe_data; + } + + virtual inline void operator() (ContextType& c, vector_type& res, matrix_type* jac = nullptr) const { + + const typename fe_var_type::fe_shape_data_type + &fe = _fe_var_data->get_fe_shape_data(); + + typename EigenVector::type + epsilon = EigenVector::type::Zero(3), + stress = EigenVector::type::Zero(3), + vec = EigenVector::type::Zero(2*fe.n_basis()); + + typename EigenMatrix::type + mat = EigenMatrix::type::Zero(3, 3), + mat1 = EigenMatrix::type::Zero(3, 2*fe.n_basis()), + mat2 = EigenMatrix::type::Zero(2*fe.n_basis(), 2*fe.n_basis()); + + MAST::FEMOperatorMatrixBase + Bxmat; + Bxmat.reinit(3, 2, fe.n_basis()); // three stress-strain components + + + for (uint_type i=0; ivalue(c, mat); + MAST::Element2D::linear_continuum_strain(*_fe_var_data, i, epsilon, Bxmat); + stress = mat * epsilon; + Bxmat.vector_mult_transpose(vec, stress); + res += fe.detJxW(i) * vec; + + if (jac) { + + Bxmat.left_multiply(mat1, mat); + Bxmat.right_multiply_transpose(mat2, mat1); + (*jac) += fe.detJxW(i) * mat2; + } + } + } protected: - const MAST::ComputeKernel *_material; - const fe_var_type *_fe_var_data; + const section_property_type *_property; + const fe_var_type *_fe_var_data; }; - +} } #endif // _mast_strain_energy_compute_kernel_h_ diff --git a/src/materials/isotropic_hookean_material.h b/src/materials/isotropic_hookean_material.h index 675ad040..ccb6d536 100644 --- a/src/materials/isotropic_hookean_material.h +++ b/src/materials/isotropic_hookean_material.h @@ -8,14 +8,14 @@ namespace MAST { -template +template class IsotropicHookeanMaterial: -public ComputeKernel { +public ComputeKernel { public: - using scalar_type = typename Traits::var_scalar_type; - using value_type = typename Traits::material_value_type; + using scalar_type = typename Traits::var_scalar_type; + using value_type = typename Traits::material_value_type; IsotropicHookeanMaterial(const std::string& nm): MAST::ComputeKernel(nm), @@ -25,8 +25,8 @@ public ComputeKernel { virtual ~IsotropicHookeanMaterial() { } - virtual inline void set_properties(const MAST::ComputeKernel& E, - const MAST::ComputeKernel& nu) { + inline void set_properties(const MAST::ComputeKernel& E, + const MAST::ComputeKernel& nu) { libmesh_assert_msg(!_E && !_nu, "Property kernels should be cleared before setting."); @@ -34,36 +34,30 @@ public ComputeKernel { _nu = ν } - virtual inline void execute() override { } - - virtual inline void execute_derivative(const std::vector& idx) override { + void inline operator(const ContextType& c, value_type& m) { libmesh_assert_msg(_E && _nu, "Property kernels should be provided before execute()."); - scalar_type - E = _E->value(), - nu = _nu->value(), + const scalar_type + E = _E->value(c), + nu = _nu->value(c), G = E/2./(1.+nu); - MAST::init_view(_matrix, 6, 6, true); - m(0, 0) = m(1, 1) = m(2, 2) = E*(1.-nu)/(1.-nu-2.*nu*nu); m(0, 1) = m(0, 2) = m(1, 0) = m(1, 2) = m(2, 0) = m(2, 1) = E*nu/(1.-nu-2.*nu*nu); m(3, 3) = m(4, 4) = m(5, 5) = E/2./(1.+nu); } - virtual inline value_type value() const override { + inline void execute_derivative(const ContextType& c, value_type& m) { libmesh_assert_msg(_E && _nu, "Property kernels should be provided before execute()."); - scalar_type - E = _E->value(), - nu = _nu->value(), + const scalar_type + E = _E->value(c), + nu = _nu->value(c), G = E/2./(1.+nu); - value_type matrix = - MAST::init_view(_matrix, 6, 6, true); m(0, 0) = m(1, 1) = m(2, 2) = E*(1.-nu)/(1.-nu-2.*nu*nu); @@ -71,8 +65,6 @@ public ComputeKernel { m(3, 3) = m(4, 4) = m(5, 5) = E/2./(1.+nu); } - virtual inline value_type value_derivative() const override { return _matrix;} - protected: const MAST::ComputeKernel *_E; diff --git a/src/mesh/fe.h b/src/mesh/fe.h index 10c2716b..e3af3993 100644 --- a/src/mesh/fe.h +++ b/src/mesh/fe.h @@ -8,53 +8,43 @@ namespace MAST { -template -struct EigenFEShapeDataViewTraits { - - using xyz_view_type = Eigen::Matrix; - using detJ_view_type = Eigen::Matrix; - using detJxW_view_type = Eigen::Matrix; - using dx_dxi_view_type = Eigen::Matrix; - using dxi_dx_view_type = Eigen::Matrix; - using dphi_dx_view_type = Eigen::Matrix; -}; - - - -template -struct EigenFESolDataViewTraits { - - using coefficient_view_type = Eigen::Matrix; - using u_view_type = Eigen::Matrix; - using du_dx_view_type = Eigen::Matrix; -}; - - - -template +template class FEBasis: ComputeKernelBase { public: using basis_scalar_type = ScalarType; - FEBasis(const std::string& nm): - MAST::ComputeKernelBase(nm), + FEBasis(const std::string& nm, const bool executable): + MAST::ComputeKernelBase(nm, executable), _quadrature (nullptr) { } virtual ~FEBasis() {} virtual inline void set_quadrature(const MAST::Quadrature& q) { _quadrature = &q;} + virtual inline uint_type dim() const = 0; virtual inline uint_type order() const = 0; virtual inline uint_type n_basis() const = 0; + + virtual inline uint_type n_q_points() const { + + libmesh_assert_msg(_quadrature, "Quadrature not initialized."); + return _quadrature->n_points(); + } + virtual inline ScalarType qp_coord(uint_type qp, uint_type x_i) const - { return _quadrature->qp_coord(qp, x_i); } + { + libmesh_assert_msg(_quadrature, "Quadrature not initialized."); + return _quadrature->qp_coord(qp, x_i); + } + virtual inline ScalarType phi(uint_type qp, uint_type phi_i) const = 0; - virtual inline ScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + virtual inline ScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + protected: const MAST::Quadrature *_quadrature; @@ -62,16 +52,16 @@ class FEBasis: ComputeKernelBase { /*! serves as a wrapper around libMesh FEBase object to provide */ -template -class libMeshFE: public FEBasis { - +template +class libMeshFE: public FEBasis { + public: - using scalar_type = Real; - using basis_scalar_type = typename MAST::FEBasis::basis_scalar_type; + using scalar_type = ScalarType; + using basis_scalar_type = typename MAST::FEBasis::basis_scalar_type; libMeshFE(const std::string& nm): - MAST::FEBasis(nm), + MAST::FEBasis(nm, false), _compute_dphi_dxi (false), _fe (nullptr) { } @@ -80,27 +70,52 @@ class libMeshFE: public FEBasis { inline void set_compute_dphi_dxi(bool f) { _compute_dphi_dxi = f;} virtual inline void reinit(const libMesh::FEBase& fe) { + _fe = &fe; const libMesh::FEMap& m = fe.get_fe_map(); - this->_dphi_dxi = - {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; + + if (_compute_dphi_dxi) + this->_dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; + else + _dphi_dxi.clear(); } - - virtual inline void execute(ContextType& c) override { } - virtual inline uint_type dim() const override { return _fe->get_dim();} + virtual inline uint_type dim() const override { + + libmesh_assert_msg(_fe, "FE not initialized."); + + return _fe->get_dim(); + } - virtual inline uint_type order() const override {return _fe->get_order();} + virtual inline uint_type order() const override { + + libmesh_assert_msg(_fe, "FE not initialized."); + + return _fe->get_order(); + } - virtual inline uint_type n_basis() const override { return _fe->n_shape_functions();} + virtual inline uint_type n_basis() const override { + + libmesh_assert_msg(_fe, "FE not initialized."); + + return _fe->n_shape_functions(); + } virtual inline basis_scalar_type phi(uint_type qp, uint_type phi_i) const override - {return _fe->get_phi()[phi_i][qp];} + { + libmesh_assert_msg(_fe, "FE not initialized."); + + return _fe->get_phi()[phi_i][qp]; + } virtual inline basis_scalar_type dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const override - { return (*_dphi_dxi[xi_i])[phi_i][qp];} - + { + libmesh_assert_msg(_compute_dphi_dxi, "FE not initialized with derivatives."); + + return (*_dphi_dxi[xi_i])[phi_i][qp]; + } + protected: bool _compute_dphi_dxi; diff --git a/src/mesh/fe_basis_derived_data.h b/src/mesh/fe_basis_derived_data.h index c4c2884d..7f394015 100644 --- a/src/mesh/fe_basis_derived_data.h +++ b/src/mesh/fe_basis_derived_data.h @@ -9,93 +9,87 @@ namespace MAST { /*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, this is any non-Lagrange basis.*/ -template -class FEBasisDerivedData: public MAST::FEShapeDataBase { - -public: + +template +class FEBasisDerivedData: +public MAST::FEShapeDataBase { - using fe_geom_derived_type = MAST::FEGeometryBasisDerivedData; +public: FEBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm), + MAST::FEShapeDataBase(nm, false), _fe_geom(nullptr) {} virtual ~FEBasisDerivedData() {} - virtual void inline set_fe_geom_derived_data(const fe_geom_derived_type& fe_geom) - { _fe_geom = &fe_geom;} - - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _fe_geom->xyz(qp, x_i);} - - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _fe_geom->detJ(qp);} - - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _fe_geom->detJxW(qp);} - - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dx_dxi(qp, x_i, xi_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dxi_dx(qp, x_i, xi_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, phi_i, x_i);} - -protected: - - const fe_geom_derived_type *_fe_geom; - - typename ViewTraits::dphi_dx_view_type _dphi_dx; }; /*! This provides the derivative of shape functions when the FE basis is different from that used for interpolation of geometry. Typically, this is any non-Lagrange basis.*/ -template -class FEBasisDerivedData>: -public MAST::FEShapeDataBase { +template +class FEBasisDerivedData: +public MAST::FEShapeDataBase { public: - using fe_derived_data_type = MAST::FEGeometryBasisDerivedData>; + using fe_derived_data_type = MAST::FEGeometryBasisDerivedData; FEBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm), + MAST::FEShapeDataBase(nm), _fe_geom (nullptr) {} virtual ~FEBasisDerivedData() {} - + virtual void inline set_fe_geom_derived_data(const fe_derived_data_type& fe_geom) - { _fe_geom = &fe_geom;} - + { + libmesh_assert_msg(!_fe_geom, "FE object already initialized"); + + _fe_geom = &fe_geom; + } + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _fe_geom->xyz(qp, x_i);} + { + libmesh_assert_msg(_fe_geom, "FE Geometry Basis not initialized."); + + return _fe_geom->xyz(qp, x_i); + } virtual inline NodalScalarType detJ(uint_type qp) const override - { return _fe_geom->detJ(qp);} + { + libmesh_assert_msg(_fe_geom, "FE Geometry Basis not initialized."); + + return _fe_geom->detJ(qp); + } virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _fe_geom->detJxW(qp);} - + { + libmesh_assert_msg(_fe_geom, "FE Geometry Basis not initialized."); + + return _fe_geom->detJxW(qp); + } + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dx_dxi(qp, x_i, xi_i);} - + { + libmesh_assert_msg(_fe_geom, "FE Geometry Basis not initialized."); + + return _fe_geom->dx_dxi(qp, x_i, xi_i); + } + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _fe_geom->dxi_dx(qp, x_i, xi_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} - + { + libmesh_assert_msg(_fe_geom, "FE Geometry Basis not initialized."); + + return _fe_geom->dxi_dx(qp, x_i, xi_i); + } + + protected: const fe_derived_data_type& _fe_geom; - - typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; }; } diff --git a/src/mesh/fe_compute_kernel.h b/src/mesh/fe_compute_kernel.h index 95becb0f..1bb5c2bf 100644 --- a/src/mesh/fe_compute_kernel.h +++ b/src/mesh/fe_compute_kernel.h @@ -7,8 +7,8 @@ namespace MAST { -template -class FEComputeKernel: public MAST::ComputeKernel { +template +class FEComputeKernel: public MAST::ComputeKernel { public: @@ -17,7 +17,7 @@ class FEComputeKernel: public MAST::ComputeKernel { using sol_scalar_type = typename Traits::sol_scalar_type; FEComputeKernel(const std::string& nm): - MAST::ComputeKernel(nm) { } + MAST::ComputeKernel(nm, false) { } virtual ~FEComputeKernel() { } }; diff --git a/src/mesh/fe_geometry_basis_derived_data.h b/src/mesh/fe_geometry_basis_derived_data.h index 67c49705..360735d2 100644 --- a/src/mesh/fe_geometry_basis_derived_data.h +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -10,42 +10,18 @@ namespace MAST { /*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate nodal locations. Typically, Lagrange shape functoins are used for this purpose */ -template +template class FEGeometryBasisDerivedData: -public MAST::FEShapeDataBase { +public MAST::FEShapeDataBase { public: FEGeometryBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm) {} + MAST::FEShapeDataBase(nm) {} virtual ~FEGeometryBasisDerivedData() {} - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _xyz(qp, x_i);} - - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _detJ(qp);} - - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _detJxW(qp);} - - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dx_dxi(qp, x_i, xi_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dxi_dx(qp, x_i, xi_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, phi_i, x_i);} - protected: - typename ViewTraits::xyz_view_type _xyz; - typename ViewTraits::detJ_view_type _detJ; - typename ViewTraits::detJxW_view_type _detJxW; - typename ViewTraits::dx_dxi_view_type _dx_dxi; - typename ViewTraits::dxi_dx_view_type _dxi_dx; - typename ViewTraits::dphi_dx_view_type _dphi_dx; }; @@ -53,42 +29,36 @@ public MAST::FEShapeDataBase { /*! This provides the derivative of shape functions when the FE basis also forms the basis for geometry interpolation used to interpolate nodal locations. Typically, Lagrange shape functoins are used for this purpose */ template -class FEGeometryBasisDerivedData, ContextType>: -public MAST::FEShapeDataBase { +class FEGeometryBasisDerivedData: +public MAST::FEShapeDataBase { public: - + FEGeometryBasisDerivedData(const std::string& nm): - MAST::FEShapeDataBase(nm) {} + MAST::FEShapeDataBase(nm) {} virtual ~FEGeometryBasisDerivedData() {} - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const override - { return _xyz(qp, x_i);} + virtual inline void reinit(const ContextType& c) override { + + uint_type + nq = this->n_q_points(); + const int_type + d = this->spatial_dim(); + + this->_xyz = EigenMatrix::type::Zero(nq, d); + this->_detJ = EigenVector::type::Zero(nq); + this->_detJxW = EigenVector::type::Zero(nq); + this->_dx_dxi = EigenMatrix::type::Zero(nq, d*d); + this->_dxi_dx = EigenMatrix::type::Zero(nq, d*d); + this->_dphi_dx = EigenMatrix::type::Zero(nq, d*d*this->n_basis()); + } + + virtual inline void reinit_for_side(const ContextType& c, uint_type s) override { } + - virtual inline NodalScalarType detJ(uint_type qp) const override - { return _detJ(qp);} - - virtual inline NodalScalarType detJxW(uint_type qp) const override - { return _detJxW(qp);} - - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i);} - - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const override - { return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i);} - - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const override - { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} - protected: - typename EigenFEShapeDataViewTraits::xyz_view_type _xyz; - typename EigenFEShapeDataViewTraits::detJ_view_type _detJ; - typename EigenFEShapeDataViewTraits::detJxW_view_type _detJxW; - typename EigenFEShapeDataViewTraits::dx_dxi_view_type _dx_dxi; - typename EigenFEShapeDataViewTraits::dxi_dx_view_type _dxi_dx; - typename EigenFEShapeDataViewTraits::dphi_dx_view_type _dphi_dx; }; } diff --git a/src/mesh/fe_shape_data_base.h b/src/mesh/fe_shape_data_base.h index 244ca124..68846489 100644 --- a/src/mesh/fe_shape_data_base.h +++ b/src/mesh/fe_shape_data_base.h @@ -7,13 +7,30 @@ namespace MAST { -template +template class FEShapeDataBase: public MAST::ComputeKernelBase { + +public: + + FEShapeDataBase(const std::string& nm): + MAST::ComputeKernelBase(nm, false) {} + + virtual ~FEShapeDataBase() {} + +}; + + + +template +class FEShapeDataBase: +public MAST::ComputeKernelBase { public: + + using dphi_dx_type = typename Eigen::Map::type>; FEShapeDataBase(const std::string& nm): - MAST::ComputeKernelBase(nm), + MAST::ComputeKernelBase(nm, false), _compute_xyz (false), _compute_Jac (false), _compute_detJ (false), @@ -22,7 +39,8 @@ class FEShapeDataBase: public MAST::ComputeKernelBase { _compute_normal (false), _spatial_dim (0), _fe_basis (nullptr) - {} + { } + virtual ~FEShapeDataBase() {} inline void set_compute_xyz(bool f) { _compute_xyz = f;} inline void set_compute_Jac(bool f) { _compute_Jac = f;} @@ -30,24 +48,79 @@ class FEShapeDataBase: public MAST::ComputeKernelBase { inline void set_compute_detJxW(bool f) { _compute_JxW = f;} inline void set_compute_dphi_dx(bool f) { _compute_dphi_dx = f;} inline void set_compute_normal(bool f) { _compute_normal = f;} - inline void set_fe_basis(MAST::FEBasis& basis) - { _fe_basis = &basis;} - virtual inline void execute(ContextType& c) override { } - //virtual inline void reinit_for_side(MAST::FEBasis& basis, uint_type s) = 0; - virtual inline uint_type ref_dim() const { return _fe_basis->dim();} - virtual inline uint_type spatial_dim() const { return _spatial_dim;} - virtual inline uint_type order() const { return _fe_basis->order();} - virtual inline uint_type n_basis() const { return _fe_basis->n_basis();} + inline void set_fe_basis(MAST::FEBasis& basis) + { + libmesh_assert_msg(!_fe_basis, "FE Basis already initialized."); + + _fe_basis = &basis; + } + + virtual inline void reinit(const ContextType& c) = 0; + virtual inline void reinit_for_side(const ContextType& c, uint_type s) = 0; + + virtual inline uint_type ref_dim() const { + + libmesh_assert_msg(_fe_basis, "FE Basis not initialized."); + return _fe_basis->dim(); + } + + virtual inline uint_type spatial_dim() const { + + libmesh_assert_msg(_fe_basis, "FE Basis not initialized."); + return _spatial_dim; + } + + virtual inline uint_type order() const { + + libmesh_assert_msg(_fe_basis, "FE Basis not initialized."); + return _fe_basis->order(); + } + + virtual inline uint_type n_basis() const { + + libmesh_assert_msg(_fe_basis, "FE Basis not initialized."); + return _fe_basis->n_basis(); + } + virtual inline BasisScalarType phi(uint_type qp, uint_type phi_i) const - { return _fe_basis->phi(qp, phi_i);} + { + libmesh_assert_msg(_fe_basis, "FE Basis not initialized."); + return _fe_basis->phi(qp, phi_i); + } + virtual inline BasisScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const - { return _fe_basis->dphi_dxi(qp, phi_i, xi_i);} - virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const = 0; - virtual inline NodalScalarType detJ(uint_type qp) const = 0; - virtual inline NodalScalarType detJxW(uint_type qp) const = 0; - virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; - virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const = 0; - virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; + { + + libmesh_assert_msg(_fe_basis, "FE Basis not initialized."); + + return _fe_basis->dphi_dxi(qp, phi_i, xi_i); + } + + virtual inline uint_type n_q_points() const { return _fe_basis->n_q_points();} + + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const + { return _xyz(qp, x_i);} + + virtual inline NodalScalarType detJ(uint_type qp) const + { return _detJ(qp);} + + virtual inline NodalScalarType detJxW(uint_type qp) const + { return _detJxW(qp);} + + virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const + { return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i);} + + virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const + { return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i);} + + inline const dphi_dx_type + dphi_dx(uint_type qp, uint_type x_i) const + { + return dphi_dx_type(_dphi_dx.row(qp).segment(x_i*this->spatial_dim(), this->n_basis()).data(), this->n_basis()); + } + + virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const + { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} protected: @@ -59,7 +132,14 @@ class FEShapeDataBase: public MAST::ComputeKernelBase { bool _compute_normal; uint_type _spatial_dim; - MAST::FEBasis *_fe_basis; + MAST::FEBasis *_fe_basis; + + typename EigenMatrix::type _xyz; + typename EigenVector::type _detJ; + typename EigenVector::type _detJxW; + typename EigenMatrix::type _dx_dxi; + typename EigenMatrix::type _dxi_dx; + typename EigenMatrix::type _dphi_dx; }; } diff --git a/src/mesh/fe_var_data.h b/src/mesh/fe_var_data.h index db37dbc6..7ffa11fb 100644 --- a/src/mesh/fe_var_data.h +++ b/src/mesh/fe_var_data.h @@ -9,72 +9,62 @@ namespace MAST { /*! provides access to the element solution vector through a memory view. */ -template +template class FEVarData: public MAST::ComputeKernelBase { public: - using var_scalar_type = typename MAST::DeducedScalarType::type; - using coefficient_view_type = typename ViewTraits::coefficient_view_type; - using fe_shape_data_type = typename MAST::FEShapeDataBase; - - FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm) {} + FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm, false) {} virtual ~FEVarData() {} - virtual inline void execute(ContextType& c) override {} - inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - inline void - set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} - inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = coeffs;} - inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} - inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) { return _du_dx(qp, comp, x_i);} -protected: - - bool *_compute_du_dx; - const fe_shape_data_type *_fe; - const typename ViewTraits::coefficient_view_type _coeffs; - typename ViewTraits::u_view_type _u; - typename ViewTraits::du_dx_view_type _du_dx; }; template -class FEVarData, ContextType>: +class FEVarData: public MAST::ComputeKernelBase { public: using var_scalar_type = typename MAST::DeducedScalarType::type; - using coefficient_view_type = typename EigenFESolDataViewTraits::coefficient_view_type; - using u_view_type = typename MAST::EigenFESolDataViewTraits::u_view_type; - using du_dx_view_type = typename MAST::EigenFESolDataViewTraits::du_dx_view_type; - using fe_shape_data_type = typename MAST::FEShapeDataBase; + using vector_type = typename EigenTraits::vector_type; + using matrix_type = typename EigenTraits::vector_type; + using fe_shape_data_type = typename MAST::FEShapeDataBase; FEVarData(const std::string& nm): - MAST::ComputeKernelBase(nm), + MAST::ComputeKernelBase(nm, false), _compute_du_dx (nullptr), _fe (nullptr), _coeffs (nullptr) {} virtual ~FEVarData() {} - virtual inline void execute(ContextType& c) override {} - inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - inline void set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} - inline void set_fe_coefficient_view(const coefficient_view_type& coeffs) { _coeffs = &coeffs;} - inline var_scalar_type u(uint_type qp, uint_type comp) { return _u(qp, comp);} - inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) + virtual inline void init() { + + } + + virtual inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + virtual inline void set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} + virtual inline const fe_shape_data_type& get_fe_shape_data() const { + + libmesh_assert_msg(_fe, "FE pointer not initialized."); + return *_fe; + } + virtual inline void set_fe_coefficient_view(const vector_type& coeffs) { _coeffs = &coeffs;} + virtual inline uint_type n_q_points() const { return _u.size();} + virtual inline var_scalar_type u(uint_type qp, uint_type comp) const { return _u(qp, comp);} + virtual inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) const { return _du_dx(qp, _fe->spatial_dim()*comp + x_i);} protected: bool _compute_du_dx; const fe_shape_data_type *_fe; - const coefficient_view_type *_coeffs; - u_view_type _u; - du_dx_view_type _du_dx; + const vector_type *_coeffs; + matrix_type _u; + matrix_type _du_dx; }; } diff --git a/src/mesh/quadrature.h b/src/mesh/quadrature.h index 6b837f73..b0978db2 100644 --- a/src/mesh/quadrature.h +++ b/src/mesh/quadrature.h @@ -16,8 +16,8 @@ class Quadrature: public MAST::ComputeKernelBase { using scalar_type = ScalarType; - Quadrature(const std::string& nm): - MAST::ComputeKernelBase(nm) {} + Quadrature(const std::string& nm, bool executable): + MAST::ComputeKernelBase(nm, executable) {} virtual ~Quadrature() {} virtual inline uint_type dim() const = 0; virtual inline uint_type order() const = 0; @@ -30,7 +30,7 @@ class Quadrature: public MAST::ComputeKernelBase { /*! This provides a */ -template +template class MappedQuadrature: public MAST::Quadrature { public: @@ -41,9 +41,9 @@ class MappedQuadrature: public MAST::Quadrature { MAST::Quadrature(nm), _dim(0) {} virtual ~MappedQuadrature() {} virtual void set_data(const uint_type dim, - const typename ViewTraits::quadrature_point_type points, - const typename ViewTraits::quadrature_point_type weights) { - _points = points; + const typename Eigen::Matrix& points, + const typename Eigen::Matrix& weights) { + _points = points; _weights = weights; } virtual inline uint_type dim() const override { return _dim;} @@ -55,8 +55,8 @@ class MappedQuadrature: public MAST::Quadrature { protected: uint_type _dim; - typename ViewTraits::quadrature_point_type _points; - typename ViewTraits::quadrature_point_type _weights; + Eigen::Matrix _points; + Eigen::Matrix _weights; }; @@ -71,22 +71,53 @@ class libMeshQuadrature: public MAST::Quadrature { /*! the quadrature object \p q is expected to be initialized outside of this class. */ libMeshQuadrature(const std::string& nm): - MAST::Quadrature(nm), + MAST::Quadrature(nm, false), _q (nullptr) {} + virtual ~libMeshQuadrature() {} - virtual inline void execute(ContextType& c) override { } - const libMesh::QBase& get_libmesh_object() const { return *_q;} - virtual inline uint_type dim() const override { return _q->get_dim();} - virtual inline uint_type order() const override { return _q->get_order();} - virtual inline uint_type n_points() const override { return _q->n_points();} - virtual inline scalar_type qp_coord(uint_type qp, uint_type xi_i) const override - { return _q->get_points()[qp](xi_i);} - virtual inline scalar_type weight(uint_type qp) const override { return _q->w(qp);} + + virtual inline void set_quadrature(const libMesh::QBase& q) { _q = &q; } + + const libMesh::QBase& get_libmesh_object() const { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return *_q; + } + + virtual inline uint_type dim() const override { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return _q->get_dim(); + } + + virtual inline uint_type order() const override { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return _q->get_order(); + } + + virtual inline uint_type n_points() const override { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return _q->n_points(); + } + + virtual inline scalar_type qp_coord(uint_type qp, uint_type xi_i) const override { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return _q->get_points()[qp](xi_i); + } + + virtual inline scalar_type weight(uint_type qp) const override { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return _q->w(qp); + } protected: - libMesh::QBase* _q; + const libMesh::QBase* _q; }; } diff --git a/src/numerics/fem_operator_matrix.cpp b/src/numerics/fem_operator_matrix.cpp index 1e8377b2..95b13e0a 100644 --- a/src/numerics/fem_operator_matrix.cpp +++ b/src/numerics/fem_operator_matrix.cpp @@ -21,18 +21,5 @@ #include "numerics/fem_operator_matrix.h" -MAST::FEMOperatorMatrix::FEMOperatorMatrix(): -_n_interpolated_vars(0), -_n_discrete_vars(0), -_n_dofs_per_var(0) -{ - -} - - -MAST::FEMOperatorMatrix::~FEMOperatorMatrix() -{ - this->clear(); -} diff --git a/src/numerics/fem_operator_matrix.h b/src/numerics/fem_operator_matrix.h index 4d8bc56c..7bd0ad04 100644 --- a/src/numerics/fem_operator_matrix.h +++ b/src/numerics/fem_operator_matrix.h @@ -31,148 +31,175 @@ namespace MAST { + +template +class FEMOperatorMatrixBase +{ +public: + FEMOperatorMatrixBase(); - class FEMOperatorMatrix - { - public: - FEMOperatorMatrix(); - - - virtual ~FEMOperatorMatrix(); - - - /*! - * clears the data structures - */ - void clear(); - - - unsigned int m() const {return _n_interpolated_vars;} - - unsigned int n() const {return _n_discrete_vars*_n_dofs_per_var;} - - void print(std::ostream& o); - - /*! - * this initializes the operator for number of rows and variables, assuming - * that all variables has the same number of dofs. This is typically the case - * for structural strain operator matrices. Note that when this method is used - * the user must set the matrix entries by calling set_shape_functions - */ - void reinit(unsigned int n_interpolated_vars, - unsigned int n_discrete_vars, - unsigned int n_discrete_dofs_per_var); - - /*! - * sets the shape function values for the block corresponding to - * \p interpolated_var and \p discrete_var. This means that the row - * \p interpolated_var, the value in columns - * \p discrete_vars*n_discrete_dofs_per_var - (discrete_vars+1)*n_discrete_dofs_per_var-1) - * will be set equal to \p shape_func . - */ - void set_shape_function(unsigned int interpolated_var, - unsigned int discrete_var, - const RealVectorX& shape_func); - - /*! - * this initializes all variables to use the same interpolation function. - * It is assumed that the number of discrete vars is same as the number of - * interpolated vars. This is typically the case for fluid elements and - * for structural element inertial matrix calculations - */ - void reinit(unsigned int n_interpolated_vars, - const RealVectorX& shape_func); - - /*! - * res = [this] * v - */ - template - void vector_mult(T& res, const T& v) const; - - - /*! - * res = v^T * [this] - */ - template - void vector_mult_transpose(T& res, const T& v) const; - - - /*! - * [R] = [this] * [M] - */ - template - void right_multiply(T& r, const T& m) const; - - - /*! - * [R] = [this]^T * [M] - */ - template - void right_multiply_transpose(T& r, const T& m) const; - - - /*! - * [R] = [this]^T * [M] - */ - template - void right_multiply_transpose(T& r, const MAST::FEMOperatorMatrix& m) const; - - - /*! - * [R] = [M] * [this] - */ - template - void left_multiply(T& r, const T& m) const; - - - /*! - * [R] = [M] * [this]^T - */ - template - void left_multiply_transpose(T& r, const T& m) const; - - - protected: - - /*! - * number of rows of the operator - */ - unsigned int _n_interpolated_vars; - - /*! - * number of discrete variables in the system - */ - unsigned int _n_discrete_vars; - - /*! - * number of dofs for each variable - */ - unsigned int _n_dofs_per_var; - - /*! - * stores the shape function values that defines the coupling - * of i_th interpolated var and j_th discrete var. Stored in - * column major format. nullptr, if values are zero, otherwise the - * value is set in the vector. - */ - std::vector _var_shape_functions; - }; + + virtual ~FEMOperatorMatrixBase(); + + + /*! + * clears the data structures + */ + void clear(); + + + unsigned int m() const {return _n_interpolated_vars;} + + unsigned int n() const {return _n_discrete_vars*_n_dofs_per_var;} + + void print(std::ostream& o); + + /*! + * this initializes the operator for number of rows and variables, assuming + * that all variables has the same number of dofs. This is typically the case + * for structural strain operator matrices. Note that when this method is used + * the user must set the matrix entries by calling set_shape_functions + */ + void reinit(unsigned int n_interpolated_vars, + unsigned int n_discrete_vars, + unsigned int n_discrete_dofs_per_var); + + /*! + * sets the shape function values for the block corresponding to + * \p interpolated_var and \p discrete_var. This means that the row + * \p interpolated_var, the value in columns + * \p discrete_vars*n_discrete_dofs_per_var - (discrete_vars+1)*n_discrete_dofs_per_var-1) + * will be set equal to \p shape_func . + */ + template + void set_shape_function(unsigned int interpolated_var, + unsigned int discrete_var, + const VecType& shape_func); + + /*! + * this initializes all variables to use the same interpolation function. + * It is assumed that the number of discrete vars is same as the number of + * interpolated vars. This is typically the case for fluid elements and + * for structural element inertial matrix calculations + */ + template + void reinit(unsigned int n_interpolated_vars, + const VecType& shape_func); + + /*! + * res = [this] * v + */ + template + void vector_mult(T& res, const T& v) const; + + + /*! + * res = v^T * [this] + */ + template + void vector_mult_transpose(T& res, const T& v) const; + + + /*! + * [R] = [this] * [M] + */ + template + void right_multiply(T& r, const T& m) const; + + + /*! + * [R] = [this]^T * [M] + */ + template + void right_multiply_transpose(T& r, const T& m) const; + + + /*! + * [R] = [this]^T * [M] + */ + template + void right_multiply_transpose(T& r, + const MAST::FEMOperatorMatrixBase& m) const; + + + /*! + * [R] = [M] * [this] + */ + template + void left_multiply(T& r, const T& m) const; + + + /*! + * [R] = [M] * [this]^T + */ + template + void left_multiply_transpose(T& r, const T& m) const; + + +protected: + + /*! + * number of rows of the operator + */ + unsigned int _n_interpolated_vars; + + /*! + * number of discrete variables in the system + */ + unsigned int _n_discrete_vars; + + /*! + * number of dofs for each variable + */ + unsigned int _n_dofs_per_var; + + /*! + * stores the shape function values that defines the coupling + * of i_th interpolated var and j_th discrete var. Stored in + * column major format. nullptr, if values are zero, otherwise the + * value is set in the vector. + */ + std::vector _var_shape_functions; +}; + + +class FEMOperatorMatrix: public MAST::FEMOperatorMatrixBase { +public: +FEMOperatorMatrix(): MAST::FEMOperatorMatrixBase() {} +virtual ~FEMOperatorMatrix() {} +}; + +} + +template +MAST::FEMOperatorMatrixBase::FEMOperatorMatrixBase(): +_n_interpolated_vars(0), +_n_discrete_vars(0), +_n_dofs_per_var(0) +{ } +template +MAST::FEMOperatorMatrixBase::~FEMOperatorMatrixBase() +{ + this->clear(); +} + +template inline void -MAST::FEMOperatorMatrix::print(std::ostream& o) { +MAST::FEMOperatorMatrixBase::print(std::ostream& o) { unsigned int index = 0; - + for (unsigned int i=0; i<_n_interpolated_vars; i++) {// row for (unsigned int j=0; j<_n_discrete_vars; j++) { // column index = j*_n_interpolated_vars+i; if (_var_shape_functions[index]) // check if this is non-nullptr for (unsigned int k=0; k<_n_dofs_per_var; k++) - o << std::setw(15) << (*_var_shape_functions[index])(k); + o << std::setw(15) << _var_shape_functions[index][k]; else for (unsigned int k=0; k<_n_dofs_per_var; k++) o << std::setw(15) << 0.; @@ -182,16 +209,18 @@ MAST::FEMOperatorMatrix::print(std::ostream& o) { } +template inline void -MAST::FEMOperatorMatrix::clear() { +MAST::FEMOperatorMatrixBase::clear() { _n_interpolated_vars = 0; _n_discrete_vars = 0; _n_dofs_per_var = 0; // iterate over the shape function entries and delete the non-nullptr values - std::vector::iterator it = _var_shape_functions.begin(), + typename std::vector::iterator + it = _var_shape_functions.begin(), end = _var_shape_functions.end(); for ( ; it!=end; it++) @@ -204,9 +233,10 @@ MAST::FEMOperatorMatrix::clear() { +template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: reinit(unsigned int n_interpolated_vars, unsigned int n_discrete_vars, unsigned int n_discrete_dofs_per_var) { @@ -215,19 +245,19 @@ reinit(unsigned int n_interpolated_vars, _n_interpolated_vars = n_interpolated_vars; _n_discrete_vars = n_discrete_vars; _n_dofs_per_var = n_discrete_dofs_per_var; - _var_shape_functions.resize(_n_interpolated_vars*_n_discrete_vars); - for (unsigned int i=0; i<_var_shape_functions.size(); i++) - _var_shape_functions[i] = nullptr; + _var_shape_functions.resize(_n_interpolated_vars*_n_discrete_vars, nullptr); } +template +template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: set_shape_function(unsigned int interpolated_var, unsigned int discrete_var, - const RealVectorX& shape_func) { + const VecType& shape_func) { // make sure that reinit has been called. libmesh_assert(_var_shape_functions.size()); @@ -235,48 +265,53 @@ set_shape_function(unsigned int interpolated_var, // also make sure that the specified indices are within bounds libmesh_assert(interpolated_var < _n_interpolated_vars); libmesh_assert(discrete_var < _n_discrete_vars); + libmesh_assert_equal_to(shape_func.size(), _n_dofs_per_var); - RealVectorX* vec = + ScalarType* vec = _var_shape_functions[discrete_var*_n_interpolated_vars+interpolated_var]; - if (!vec) - { - vec = new RealVectorX; - _var_shape_functions[discrete_var*_n_interpolated_vars+interpolated_var] = vec; - } + libmesh_assert_msg(!vec, "Shape function already set"); + + vec = new ScalarType[shape_func.size()]; + _var_shape_functions[discrete_var*_n_interpolated_vars+interpolated_var] = vec; - *vec = shape_func; + for (uint_type i=0; i<_n_dofs_per_var; i++) + vec[i] = shape_func(i); } +template +template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: reinit(unsigned int n_vars, - const RealVectorX& shape_func) { + const VecType& shape_func) { this->clear(); _n_interpolated_vars = n_vars; _n_discrete_vars = n_vars; _n_dofs_per_var = (unsigned int)shape_func.size(); - _var_shape_functions.resize(n_vars*n_vars); + _var_shape_functions.resize(n_vars*n_vars, nullptr); for (unsigned int i=0; i template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: vector_mult(T& res, const T& v) const { libmesh_assert_equal_to(res.size(), _n_interpolated_vars); @@ -291,15 +326,15 @@ vector_mult(T& res, const T& v) const { if (_var_shape_functions[index]) // check if this is non-nullptr for (unsigned int k=0; k<_n_dofs_per_var; k++) res(i) += - (*_var_shape_functions[index])(k) * v(j*_n_dofs_per_var+k); + _var_shape_functions[index][k] * v(j*_n_dofs_per_var+k); } } - +template template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: vector_mult_transpose(T& res, const T& v) const { libmesh_assert_equal_to(res.size(), n()); @@ -314,16 +349,17 @@ vector_mult_transpose(T& res, const T& v) const { if (_var_shape_functions[index]) // check if this is non-nullptr for (unsigned int k=0; k<_n_dofs_per_var; k++) res(j*_n_dofs_per_var+k) += - (*_var_shape_functions[index])(k) * v(i); + _var_shape_functions[index][k] * v(i); } } +template template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: right_multiply(T& r, const T& m) const { libmesh_assert_equal_to(r.rows(), _n_interpolated_vars); @@ -340,7 +376,7 @@ right_multiply(T& r, const T& m) const { for (unsigned int l=0; l template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: right_multiply_transpose(T& r, const T& m) const { libmesh_assert_equal_to(r.rows(), n()); @@ -368,18 +405,19 @@ right_multiply_transpose(T& r, const T& m) const { for (unsigned int l=0; l template inline void -MAST::FEMOperatorMatrix:: -right_multiply_transpose(T& r, const MAST::FEMOperatorMatrix& m) const { +MAST::FEMOperatorMatrixBase:: +right_multiply_transpose(T& r, const MAST::FEMOperatorMatrixBase& m) const { libmesh_assert_equal_to(r.rows(), n()); libmesh_assert_equal_to(r.cols(), m.n()); @@ -395,12 +433,13 @@ right_multiply_transpose(T& r, const MAST::FEMOperatorMatrix& m) const { index_j = j*m._n_interpolated_vars+k; if (_var_shape_functions[index_i] && m._var_shape_functions[index_j]) { // if shape function exists for both - const RealVectorX &n1 = *_var_shape_functions[index_i], - &n2 = *m._var_shape_functions[index_j]; - for (unsigned int i_n1=0; i_n1 template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: left_multiply(T& r, const T& m) const { libmesh_assert_equal_to(r.rows(), m.rows()); @@ -428,17 +468,18 @@ left_multiply(T& r, const T& m) const { for (unsigned int l=0; l template inline void -MAST::FEMOperatorMatrix:: +MAST::FEMOperatorMatrixBase:: left_multiply_transpose(T& r, const T& m) const { libmesh_assert_equal_to(r.rows(), m.rows()); @@ -455,13 +496,12 @@ left_multiply_transpose(T& r, const T& m) const { for (unsigned int l=0; l Date: Wed, 10 Jun 2020 21:00:30 -0500 Subject: [PATCH 7/9] -- Added implementations for multiple initializations for finite element basis functions. -- Added a new elem operations class to replace the nonlinear implicit assembly operations for structural analysis. --- src/elasticity/elasticity_elem_operations.h | 128 ++++++++---- ...structural_nonlinear_elem_compute_kernel.h | 195 ++++++++++++++++++ src/mesh/fe.h | 66 +++++- src/mesh/fe_geometry_basis_derived_data.h | 112 +++++++++- src/mesh/fe_shape_data_base.h | 82 +++++++- src/mesh/fe_var_data.h | 91 +++++++- src/mesh/quadrature.h | 53 ++++- 7 files changed, 655 insertions(+), 72 deletions(-) create mode 100644 src/elasticity/structural_nonlinear_elem_compute_kernel.h diff --git a/src/elasticity/elasticity_elem_operations.h b/src/elasticity/elasticity_elem_operations.h index 2d6b7a43..fa0b4794 100644 --- a/src/elasticity/elasticity_elem_operations.h +++ b/src/elasticity/elasticity_elem_operations.h @@ -48,8 +48,28 @@ struct ElasticityTraits { struct ElasticityElemOperationsContext { - ElasticityElemOperationsContext(): qp(0) {} + ElasticityElemOperationsContext(): + qp (0), + elem (nullptr), + sys (nullptr), + physics (nullptr), + sol (nullptr) {} + uint_type qp; + const libMesh::Elem* elem; + const MAST::NonlinearSystem* sys; + const MAST::PhysicsDisciplineBase* physics; + const libMesh::NumericVector* sol; + + inline Real nodal_coord(uint_type nd, uint_type d) const { + libmesh_assert_msg(elem, "Element not initialized."); + return elem->point(nd)(d); + } + + inline uint_type n_nodes() const { + libmesh_assert_msg(elem, "Element not initialized."); + return elem->n_nodes(); + } }; @@ -60,9 +80,7 @@ public MAST::ComputationBase { public: using context_type = typename Traits::context_type; - 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; + using var_scalar_type = typename Traits::var_scalar_type; using quadrature_type = typename Traits::quadrature_type; using fe_basis_type = typename Traits::fe_basis_type; using fe_derived_type = typename Traits::fe_derived_type; @@ -75,17 +93,78 @@ public MAST::ComputationBase { _sys (nullptr), _physics (nullptr), _X (nullptr), + _quadrature (nullptr), + _fe_basis (nullptr), + _fe_derived (nullptr), _fe_var_data (nullptr), _strain_energy (nullptr) { } - virtual ~ElasticityElemOperations() {} + + virtual ~ElasticityElemOperations() { this->clear();} + virtual void init(const MAST::NonlinearSystem& sys, const MAST::PhysicsDisciplineBase& physics, - const libMesh::NumericVector& X); - virtual void clear() {} - virtual void reinit(const libMesh::Elem& e) {} - virtual void compute_residual() {} - virtual void compute_jacobian() {} + const libMesh::NumericVector& sol) { + + libmesh_assert(!_initialized); + + + // setup the volume compute kernels + _sys = &sys; + _physics = &physics; + _X = / + + _quadrature = new quadrature_type("quadrature"); + _fe_basis = new fe_basis_type("fe_basis"); + _fe_basis->set_quadrature(*_quadrature); + _fe_derived = new fe_derived_type("fe_derived"); + _fe_derived->set_fe_basis(*_fe_basis); + _fe_var_data = new fe_var_type("u_vars"); + _fe_var_data->set_fe_shape_data(*_fe_derived); + _strain_energy = new MAST::Element2D::LinearContinuumStrainEnergy; + _strain_energy->set_fe_var_data(*_fe_var_data); + //_strain_energy->set_material(*_material); + + _quadrature->init_quadrature(libMesh::QGAUSS, 2); + libMesh::FEType fe_type(1, libMesh::LAGRANGE); + _fe_basis->init(_quadrature->get_libmesh_object(), fe_type); + + _initialized = true; + } + + virtual void clear() { + + if (_initialized) { + + _sys = nullptr; + _physics = nullptr; + delete _quadrature; _quadrature = nullptr; + delete _fe_basis; _fe_basis = nullptr; + delete _fe_derived; _fe_derived = nullptr; + delete _fe_var_data; _fe_var_data = nullptr; + delete _strain_energy; _strain_energy = nullptr; + } + } + + virtual void reinit(const libMesh::Elem& e) { + + libmesh_assert_msg(_initialized, "Object not initialized."); + + context_type c; + c.elem = &e; + c.sys = _sys; + c.physics = _physics; + c.sol = _X; + + _fe_basis->reinit(e); + _fe_derived->reinit(c); + _fe_var_data->init(c); + } + + virtual void compute(EigenVector& res, + EigenMatrix* jac) {} + virtual void compute_sensitivity(const MAST::FunctionBase& f, + EigenVector& res) {} virtual void compute_complex_step_jacobian() {} virtual void compute_auto_diff_jacobian() {} @@ -102,35 +181,6 @@ public MAST::ComputationBase { fe_var_type *_fe_var_data; typename MAST::Element2D::LinearContinuumStrainEnergy *_strain_energy; }; - - -template -void -MAST::ElasticityElemOperations:: -init(const MAST::NonlinearSystem& sys, - const MAST::PhysicsDisciplineBase& physics, - const libMesh::NumericVector& X) { - - libmesh_assert(!_initialized); - - - // setup the volume compute kernels - _sys = &sys; - _physics = &physics; - _X = &X; - - - _quadrature = new quadrature_type("quadrature"); - _fe_basis = new fe_basis_type("fe_basis"); - _fe_basis->set_quadrature(*_quadrature); - _fe_derived = new fe_derived_type("fe_derived"); - _fe_derived->set_fe_basis(*_fe_basis); - _fe_var_data = new fe_var_type("u_vars"); - _fe_var_data->set_fe_shape_data(*_fe_derived); - _strain_energy = new MAST::Element2D::LinearContinuumStrainEnergy; - _strain_energy->set_fe_var_data(*_fe_var_data); - //_strain_energy->set_material(*_material); -} } #endif // _mast_elasticity_elem_operations_h_ diff --git a/src/elasticity/structural_nonlinear_elem_compute_kernel.h b/src/elasticity/structural_nonlinear_elem_compute_kernel.h new file mode 100644 index 00000000..9d6e21c3 --- /dev/null +++ b/src/elasticity/structural_nonlinear_elem_compute_kernel.h @@ -0,0 +1,195 @@ +/* + * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit + * Copyright (C) 2013-2020 Manav Bhatia and MAST authors + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +#ifndef __mast__structural_nonlinear_assembly_elem_operations__ +#define __mast__structural_nonlinear_assembly_elem_operations__ + +// MAST includes +#include "base/nonlinear_implicit_assembly_elem_operations.h" +#include "elasticity/elasticity_elem_operations.h" + +namespace MAST { + + + // Forward declerations + class StructuralAssembly; + + + class StructuralNonlinearAssemblyComputations: + public MAST::NonlinearImplicitAssemblyElemOperations { + + public: + + /*! + * constructor associates this assembly object with the system + */ + StructuralNonlinearAssemblyComputations(): + MAST::NonlinearImplicitAssemblyElemOperations() { } + + + /*! + * destructor resets the association of this assembly object with + * the system + */ + virtual ~StructuralNonlinearAssemblyComputations() { + + delete _elem_ops; + } + + /*! + * sets the element solution(s) before calculations + */ + virtual void set_elem_solution(const RealVectorX& sol) { + + libmesh_assert_msg(_elem_ops, "Object not initialized"); + + _elem_ops->compute(vec, if_jac?&mat:nullptr); + } + + /*! + * performs the element calculations over \p elem, and returns + * the element vector and matrix quantities in \p mat and + * \p vec, respectively. \p if_jac tells the method to also + * assemble the Jacobian, in addition to the residual vector. + */ + virtual void elem_calculations(bool if_jac, + RealVectorX& vec, + RealMatrixX& mat) { + + libmesh_assert_msg(_elem_ops, "Object not initialized"); + + _elem_ops->compute(vec, if_jac?&mat:nullptr); + } + + + /*! + * performs the element calculations over \p elem, and returns + * the element vector quantity in \p vec. The vector quantity only + * include the \f$ [J] \{dX\} \f$ components, so the inherited classes + * must ensure that no component of constant forces (traction/body + * forces/etc.) are added to this vector. + */ + virtual void + elem_linearized_jacobian_solution_product(RealVectorX& vec) + { libmesh_assert_msg(false, "not implemented yet.");} + + + /*! + * performs the element sensitivity calculations over \p elem, + * and returns the element residual sensitivity in \p vec . + */ + virtual void elem_sensitivity_calculations(const MAST::FunctionBase& f, + RealVectorX& vec) { + + libmesh_assert_msg(_elem_ops, "Object not initialized"); + + _elem_ops->compute_sensitivity(f, vec); + } + + + /*! + * performs the element shape sensitivity calculations over \p elem, + * and returns the element residual sensitivity in \p vec . + */ + virtual void + elem_shape_sensitivity_calculations(const MAST::FunctionBase& f, + RealVectorX& vec) + { libmesh_assert_msg(false, "not implemented yet.");} + + + /*! + * performs the element topology sensitivity calculations over \p elem, + * and returns the element residual sensitivity in \p vec . + */ + virtual void + elem_topology_sensitivity_calculations(const MAST::FunctionBase& f, + RealVectorX& vec) + { libmesh_assert_msg(false, "not implemented yet.");} + + + + /*! + * performs the element topology sensitivity calculations over \p elem, + * and returns the element residual sensitivity in \p vec . + */ + virtual void + elem_topology_sensitivity_calculations(const MAST::FunctionBase& f, + const MAST::FieldFunction& vel, + RealVectorX& vec) + { libmesh_assert_msg(false, "not implemented yet.");} + + + /*! + * calculates \f$ d ([J] \{\Delta X\})/ dX \f$ over \p elem, + * and returns the matrix in \p vec . + */ + virtual void + elem_second_derivative_dot_solution_assembly(RealMatrixX& mat) + { libmesh_assert_msg(false, "not implemented yet.");} + + + /*! + * sets the structural element y-vector if 1D element is used. + */ + virtual void + set_elem_data(unsigned int dim, + const libMesh::Elem& ref_elem, + MAST::GeomElem& elem) const + { libmesh_assert_msg(false, "not implemented yet.");} + + + /*! + * initializes the object for the geometric element \p elem. This + * expects the object to be in a cleared state, so the user should + * call \p clear_elem() between successive initializations. + */ + virtual void + init(const MAST::GeomElem& elem) { + + libmesh_assert_msg(!_elem_ops, "Object already initialized"); + + _elem_ops = MAST::ElasticityElemOperations>; + _elem_ops->init(_system->system(), *_discipline); + } + + + /*! + * clears the element initialization + */ + virtual void clear_elem() { + + if (_elem_ops) { + + delete _elem_ops; + _elem_ops = nullptr; + } + + MAST::NonlinearImplicitAssemblyElemOperations::clear_elem(); + } + + + protected: + + + MAST::ElasticityElemOperations> *_elem_ops; + }; +} + + +#endif // __mast__structural_nonlinear_assembly_elem_operations__ diff --git a/src/mesh/fe.h b/src/mesh/fe.h index e3af3993..19d18a91 100644 --- a/src/mesh/fe.h +++ b/src/mesh/fe.h @@ -40,7 +40,13 @@ class FEBasis: ComputeKernelBase { libmesh_assert_msg(_quadrature, "Quadrature not initialized."); return _quadrature->qp_coord(qp, x_i); } - + + virtual inline ScalarType qp_weight(uint_type qp) const + { + libmesh_assert_msg(_quadrature, "Quadrature not initialized."); + return _quadrature->weight(qp); + } + virtual inline ScalarType phi(uint_type qp, uint_type phi_i) const = 0; virtual inline ScalarType dphi_dxi(uint_type qp, uint_type phi_i, uint_type xi_i) const = 0; @@ -63,16 +69,64 @@ class libMeshFE: public FEBasis { libMeshFE(const std::string& nm): MAST::FEBasis(nm, false), _compute_dphi_dxi (false), + _own_pointer (false), _fe (nullptr) { } - virtual ~libMeshFE() { } + + virtual ~libMeshFE() { + + if (_own_pointer) + delete _fe; + } inline void set_compute_dphi_dxi(bool f) { _compute_dphi_dxi = f;} - virtual inline void reinit(const libMesh::FEBase& fe) { + virtual inline void init(libMesh::FEBase& fe) { + + libmesh_assert_msg(!_fe, "Object already initialized"); - _fe = &fe; + _fe = &fe; + _own_pointer = false; + const libMesh::FEMap& m = fe.get_fe_map(); + if (_compute_dphi_dxi) + this->_dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; + else + _dphi_dxi.clear(); + } + + + virtual inline void init(libMesh::QBase& q, const libMesh::FEType fe_type) { + + libmesh_assert_msg(!_fe, "Object already initialized"); + + _fe = libMesh::FEBase::build(q.get_dim(), fe_type).release(); + _own_pointer = true; + _fe->attach_quadrature_rule(&q); + } + + + virtual inline void clear() { + + if (_own_pointer) + delete _fe; + + _own_pointer = false; + _fe = nullptr; + _dphi_dxi.clear(); + } + + + virtual inline void reinit(const libMesh::Elem& e) { + + libmesh_assert_msg(_fe, "Object not initialized"); + + if (_compute_dphi_dxi) + _fe->get_JxW(); + + _fe->reinit(&e); + + const libMesh::FEMap& m = _fe->get_fe_map(); if (_compute_dphi_dxi) this->_dphi_dxi = {&m.get_dphidxi_map(), &m.get_dphideta_map(), &m.get_dphidzeta_map()}; @@ -80,6 +134,7 @@ class libMeshFE: public FEBasis { _dphi_dxi.clear(); } + virtual inline uint_type dim() const override { libmesh_assert_msg(_fe, "FE not initialized."); @@ -119,7 +174,8 @@ class libMeshFE: public FEBasis { protected: bool _compute_dphi_dxi; - const libMesh::FEBase* _fe; + bool _own_pointer; + libMesh::FEBase* _fe; std::vector>*> _dphi_dxi; }; diff --git a/src/mesh/fe_geometry_basis_derived_data.h b/src/mesh/fe_geometry_basis_derived_data.h index 360735d2..105ee4e2 100644 --- a/src/mesh/fe_geometry_basis_derived_data.h +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -45,13 +45,113 @@ public MAST::FEShapeDataBasen_q_points(); const int_type d = this->spatial_dim(); + const uint_type + n_nodes = c.n_nodes(); - this->_xyz = EigenMatrix::type::Zero(nq, d); - this->_detJ = EigenVector::type::Zero(nq); - this->_detJxW = EigenVector::type::Zero(nq); - this->_dx_dxi = EigenMatrix::type::Zero(nq, d*d); - this->_dxi_dx = EigenMatrix::type::Zero(nq, d*d); - this->_dphi_dx = EigenMatrix::type::Zero(nq, d*d*this->n_basis()); + // for this class the number of basis functions should be equal to the number + // of nodes + libmesh_assert_equal_to(n_nodes, this->_fe_basis->n_basis()); + + if (this->_compute_xyz) { + + this->_node_coord = EigenMatrix::type::Zero(n_nodes, d); + this->_xyz = EigenMatrix::type::Zero(nq, d); + + // get the nodal locations + for (uint_type i=0; i_node_coord(i, j) = c.nodal_coord(i, j); + + // quadrature point coordinates + for (uint_type i=0; i_xyz(i, j) += this->_fe_basis->phi(i, k) * this->_node_coord(k, j); + } + + + if (this->_compute_Jac) { + + this->_dx_dxi = EigenMatrix::type::Zero(nq, d*d); + + // quadrature point spatial coordinate derivatives dx/dxi + for (uint_type i=0; i::type> + dxdxi(this->_dx_dxi.row(i).data(), d, d); + + for (uint_type l=0; l_fe_basis->dphi_dxi(i, l, k) * this->_node_coord(l, j); + } + } + + + if (this->_compute_detJ) { + + this->_detJ = EigenVector::type::Zero(nq); + + for (uint_type i=0; i::type> + dxdxi(this->_dx_dxi.row(i).data(), d, d); + + // determinant of dx_dxi + this->_detJ(i) = dxdxi.determinant(); + } + } + + + if (this->_compute_JxW) { + + this->_detJxW = EigenVector::type::Zero(nq); + + for (uint_type i=0; i::type> + dxdxi(this->_dx_dxi.row(i).data(), d, d); + + // determinant times weight + this->_detJxW(i) = this->_detJ(i) * this->_fe_basis->qp_weight(i); + } + } + + + if (this->_compute_Jac_inv) { + + this->_dxi_dx = EigenMatrix::type::Zero(nq, d*d); + + for (uint_type i=0; i::type> + dxdxi(this->_dx_dxi.row(i).data(), d, d), + dxidx(this->_dxi_dx.row(i).data(), d, d); + + // compute dx/dxi + dxidx = dxdxi.inverse(); + } + } + + + if (this->_compute_dphi_dx) { + + this->_dphi_dx = EigenMatrix::type::Zero(nq, d*d*this->n_basis()); + + for (uint_type i=0; i::type> + dxidx (this->_dxi_dx.row(i).data(), d, d), + dphidx(this->_dphi_dx.row(i).data(), d, n_nodes); + + for (uint_type l=0; l_fe_basis->dphi_dxi(i, l, k) * dxidx(k, j); + } + } } virtual inline void reinit_for_side(const ContextType& c, uint_type s) override { } diff --git a/src/mesh/fe_shape_data_base.h b/src/mesh/fe_shape_data_base.h index 68846489..e12192ce 100644 --- a/src/mesh/fe_shape_data_base.h +++ b/src/mesh/fe_shape_data_base.h @@ -33,6 +33,7 @@ public MAST::ComputeKernelBase { MAST::ComputeKernelBase(nm, false), _compute_xyz (false), _compute_Jac (false), + _compute_Jac_inv (false), _compute_detJ (false), _compute_JxW (false), _compute_dphi_dx (false), @@ -42,12 +43,45 @@ public MAST::ComputeKernelBase { { } virtual ~FEShapeDataBase() {} + inline void set_compute_xyz(bool f) { _compute_xyz = f;} - inline void set_compute_Jac(bool f) { _compute_Jac = f;} - inline void set_compute_detJ(bool f) { _compute_detJ = f;} - inline void set_compute_detJxW(bool f) { _compute_JxW = f;} - inline void set_compute_dphi_dx(bool f) { _compute_dphi_dx = f;} - inline void set_compute_normal(bool f) { _compute_normal = f;} + + inline void set_compute_Jac(bool f) { + + _compute_Jac = f; + if (f) this->set_compute_xyz(true); + } + + inline void set_compute_Jac_inverse(bool f) { + + _compute_Jac_inv = f; + if (f) this->set_compute_Jac(true); + } + + inline void set_compute_detJ(bool f) { + + _compute_detJ = f; + if (f) this->set_compute_Jac(true); + } + + inline void set_compute_detJxW(bool f) { + + _compute_JxW = f; + if (f) this->set_compute_Jac(true); + } + + inline void set_compute_dphi_dx(bool f) { + + _compute_dphi_dx = f; + if (f) this->set_compute_Jac_inverse(true); + } + + inline void set_compute_normal(bool f) { + + _compute_normal = f; + if (f) this->set_compute_Jac(true); + } + inline void set_fe_basis(MAST::FEBasis& basis) { libmesh_assert_msg(!_fe_basis, "FE Basis already initialized."); @@ -98,25 +132,49 @@ public MAST::ComputeKernelBase { virtual inline uint_type n_q_points() const { return _fe_basis->n_q_points();} + virtual inline NodalScalarType node_coord(uint_type nd, uint_type x_i) const + { + libmesh_assert_msg(_compute_xyz, "Nodal and QPoint locations not requested"); + return _node_coord(nd, x_i); + } + virtual inline NodalScalarType xyz(uint_type qp, uint_type x_i) const - { return _xyz(qp, x_i);} + { + libmesh_assert_msg(_compute_xyz, "Nodal and QPoint locations not requested"); + return _xyz(qp, x_i); + } virtual inline NodalScalarType detJ(uint_type qp) const - { return _detJ(qp);} + { + libmesh_assert_msg(_compute_detJ, "Jacobian computation not requested"); + return _detJ(qp); + } virtual inline NodalScalarType detJxW(uint_type qp) const - { return _detJxW(qp);} + { + libmesh_assert_msg(_compute_JxW, "JxW computation not requested"); + return _detJxW(qp); + } virtual inline NodalScalarType dx_dxi(uint_type qp, uint_type x_i, uint_type xi_i) const - { return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i);} + { + libmesh_assert_msg(_compute_Jac, "Jacobian computation not requested"); + return _dx_dxi(qp, xi_i*this->spatial_dim()+x_i); + } virtual inline NodalScalarType dxi_dx(uint_type qp, uint_type x_i, uint_type xi_i) const - { return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i);} + { + libmesh_assert_msg(_compute_Jac_inv, "Jacobian inverse computation not requested"); + return _dxi_dx(qp, xi_i*this->spatial_dim()+x_i); + } inline const dphi_dx_type dphi_dx(uint_type qp, uint_type x_i) const { - return dphi_dx_type(_dphi_dx.row(qp).segment(x_i*this->spatial_dim(), this->n_basis()).data(), this->n_basis()); + libmesh_assert_msg(_compute_dphi_dx, "Jacobian inverse computation not requested"); + + return dphi_dx_type(_dphi_dx.row(qp).segment(x_i*this->n_basis(), this->n_basis()).data(), + this->n_basis()); } virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const @@ -126,6 +184,7 @@ public MAST::ComputeKernelBase { bool _compute_xyz; bool _compute_Jac; + bool _compute_Jac_inv; bool _compute_detJ; bool _compute_JxW; bool _compute_dphi_dx; @@ -134,6 +193,7 @@ public MAST::ComputeKernelBase { MAST::FEBasis *_fe_basis; + typename EigenMatrix::type _node_coord; typename EigenMatrix::type _xyz; typename EigenVector::type _detJ; typename EigenVector::type _detJxW; diff --git a/src/mesh/fe_var_data.h b/src/mesh/fe_var_data.h index 7ffa11fb..ccbc7044 100644 --- a/src/mesh/fe_var_data.h +++ b/src/mesh/fe_var_data.h @@ -41,28 +41,109 @@ public MAST::ComputeKernelBase { {} virtual ~FEVarData() {} - virtual inline void init() { + virtual inline void init(const ContextType& c, + const std::vector& active_vars, + const vector_type* coeffs = nullptr) { + libmesh_assert_msg(!_coeffs, "Object already initialized"); + libmesh_assert_equal_to(_fe->spatial_dim(), c.elem->dim()); + + // first get the coefficients. If the vector coeffs is provided then + // use that, else extract it from the context. + if (coeffs) { + _coeffs = coeffs; + _coeff_vec.setZero(); + } + else { + + std::vector dofs; + c.sys->get_dof_map().dof_indices(c.elem, dofs); + + uint_type + n_dofs = dofs.size(); + + _coeff_vec = vector_type::Zero(n_dofs); + + for (uint_type i=0; iel(dofs[i]); + + _coeffs = &_coeff_vec; + } + + // number of avtive variables should be less than or equal to the number of + // variables in the system + libmesh_assert_less_equal(active_vars.size(), c.sys->n_vars()); + + _active_vars = active_vars; + + uint_type + d = _fe->spatial_dim(), + n_qp = _fe->n_q_points(), + n_basis = _fe->n_basis(), + n_comp = active_vars.size(); + + libmesh_assert_equal_to(_coeffs->size(), n_comp*n_basis); + + _u = vector_type::Zero(n_qp, n_comp); + + // now, initialize the solution value and derivatives. + for (uint_type i=0; iphi(i, j) * (*_coeffs)(j*n_comp+k); + + if (_compute_du_dx) { + + _du_dx = vector_type::Zero(n_qp, n_comp*d); + + for (uint_type i=0; idphi_dx(i, k, l) * (*_coeffs)(j*n_comp+k); + + } } virtual inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + virtual inline void set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} + virtual inline const fe_shape_data_type& get_fe_shape_data() const { libmesh_assert_msg(_fe, "FE pointer not initialized."); return *_fe; } - virtual inline void set_fe_coefficient_view(const vector_type& coeffs) { _coeffs = &coeffs;} - virtual inline uint_type n_q_points() const { return _u.size();} - virtual inline var_scalar_type u(uint_type qp, uint_type comp) const { return _u(qp, comp);} + + virtual inline uint_type n_q_points() const { + + libmesh_assert_msg(_fe, "FE pointer not initialized"); + return _fe->n_q_points(); + } + + virtual inline var_scalar_type u(uint_type qp, uint_type comp) const { + + libmesh_assert_msg(_coeffs, "Object not initialized"); + libmesh_assert_less_equal(comp, _active_vars.size()); + + return _u(qp, comp); + } + virtual inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) const - { return _du_dx(qp, _fe->spatial_dim()*comp + x_i);} + { + libmesh_assert_msg(_coeffs, "Object not initialized"); + libmesh_assert_less_equal(comp, _active_vars.size()); + + return _du_dx(qp, _fe->spatial_dim()*comp + x_i); + } protected: bool _compute_du_dx; const fe_shape_data_type *_fe; const vector_type *_coeffs; + std::vector _active_vars; + vector_type _coeff_vec; matrix_type _u; matrix_type _du_dx; }; diff --git a/src/mesh/quadrature.h b/src/mesh/quadrature.h index b0978db2..0ca62eb2 100644 --- a/src/mesh/quadrature.h +++ b/src/mesh/quadrature.h @@ -6,6 +6,10 @@ // MAST includes #include "base/compute_kernel.h" +// libMesh includes +#include "libmesh/quadrature.h" +#include "libmesh/enum_quadrature_type.h" + namespace MAST { @@ -72,13 +76,49 @@ class libMeshQuadrature: public MAST::Quadrature { /*! the quadrature object \p q is expected to be initialized outside of this class. */ libMeshQuadrature(const std::string& nm): MAST::Quadrature(nm, false), - _q (nullptr) - {} + _own_pointer (false), + _q (nullptr) + { } - virtual ~libMeshQuadrature() {} + virtual ~libMeshQuadrature() { + + if (_q && _own_pointer) + delete _q; + } - virtual inline void set_quadrature(const libMesh::QBase& q) { _q = &q; } + virtual inline void set_quadrature(libMesh::QBase& q) { + + libmesh_assert_msg(!_q, "Quadrature already initialized."); + + _q = &q; + _own_pointer = false; + } + + virtual inline void init_quadrature(const libMesh::QuadratureType qt, + const uint_type dim, + const libMesh::Order order=libMesh::INVALID_ORDER) + { + + libmesh_assert_msg(!_q, "Quadrature already initialized."); + + _q = libMesh::QBase::build (qt, dim, order).release(); + _own_pointer = true; + } + virtual inline void clear() { + + if (_q && _own_pointer) + delete _q; + + _q = nullptr; + } + + libMesh::QBase& get_libmesh_object() { + + libmesh_assert_msg(_q, "Quadrature object not initialized"); + return *_q; + } + const libMesh::QBase& get_libmesh_object() const { libmesh_assert_msg(_q, "Quadrature object not initialized"); @@ -116,8 +156,9 @@ class libMeshQuadrature: public MAST::Quadrature { } protected: - - const libMesh::QBase* _q; + + bool _own_pointer; + libMesh::QBase* _q; }; } From ed97e8df995c13679dd9f19d4e412a551474b549 Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Thu, 11 Jun 2020 14:43:20 -0500 Subject: [PATCH 8/9] -- Several bug-fixes for templated compute kernel implementation of elements and elasticity operation. --- src/elasticity/elasticity_elem_operations.h | 99 ++++++++++++------- src/elasticity/strain_energy_compute_kernel.h | 9 ++ src/mesh/fe.h | 8 +- src/mesh/fe_geometry_basis_derived_data.h | 4 +- src/mesh/fe_shape_data_base.h | 8 +- src/mesh/fe_var_data.h | 60 ++++++----- src/numerics/fem_operator_matrix.h | 9 +- 7 files changed, 127 insertions(+), 70 deletions(-) diff --git a/src/elasticity/elasticity_elem_operations.h b/src/elasticity/elasticity_elem_operations.h index fa0b4794..18aaa9b0 100644 --- a/src/elasticity/elasticity_elem_operations.h +++ b/src/elasticity/elasticity_elem_operations.h @@ -85,26 +85,29 @@ public MAST::ComputationBase { using fe_basis_type = typename Traits::fe_basis_type; using fe_derived_type = typename Traits::fe_derived_type; using fe_var_type = typename Traits::fe_var_type; + using vector_type = typename Traits::view_traits::template vector_type; + using matrix_type = typename Traits::view_traits::template matrix_type; using section_property_type = typename Traits::section_property_type; ElasticityElemOperations(): MAST::ComputationBase(), - _initialized (false), - _sys (nullptr), - _physics (nullptr), - _X (nullptr), - _quadrature (nullptr), - _fe_basis (nullptr), - _fe_derived (nullptr), - _fe_var_data (nullptr), - _strain_energy (nullptr) + _initialized (false), + _sys (nullptr), + _physics (nullptr), + _X (nullptr), + _quadrature (nullptr), + _fe_basis (nullptr), + _fe_derived (nullptr), + _fe_var_data (nullptr), + _strain_energy (nullptr), + _section_property (nullptr) { } virtual ~ElasticityElemOperations() { this->clear();} - virtual void init(const MAST::NonlinearSystem& sys, - const MAST::PhysicsDisciplineBase& physics, - const libMesh::NumericVector& sol) { + virtual inline void init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& sol) { libmesh_assert(!_initialized); @@ -114,25 +117,27 @@ public MAST::ComputationBase { _physics = &physics; _X = / - _quadrature = new quadrature_type("quadrature"); - _fe_basis = new fe_basis_type("fe_basis"); + _quadrature = new quadrature_type("quadrature"); + _fe_basis = new fe_basis_type("fe_basis"); _fe_basis->set_quadrature(*_quadrature); - _fe_derived = new fe_derived_type("fe_derived"); + _fe_derived = new fe_derived_type("fe_derived"); _fe_derived->set_fe_basis(*_fe_basis); - _fe_var_data = new fe_var_type("u_vars"); + _fe_var_data = new fe_var_type("u_vars"); _fe_var_data->set_fe_shape_data(*_fe_derived); - _strain_energy = new MAST::Element2D::LinearContinuumStrainEnergy; + _strain_energy = new MAST::Element2D::LinearContinuumStrainEnergy; _strain_energy->set_fe_var_data(*_fe_var_data); - //_strain_energy->set_material(*_material); + _section_property = new section_property_type; + _strain_energy->set_section_property(*_section_property); - _quadrature->init_quadrature(libMesh::QGAUSS, 2); + _quadrature->init_quadrature(libMesh::QGAUSS, 2, libMesh::SECOND); libMesh::FEType fe_type(1, libMesh::LAGRANGE); + _fe_basis->set_compute_dphi_dxi(true); _fe_basis->init(_quadrature->get_libmesh_object(), fe_type); _initialized = true; } - virtual void clear() { + virtual inline void clear() { if (_initialized) { @@ -146,27 +151,51 @@ public MAST::ComputationBase { } } - virtual void reinit(const libMesh::Elem& e) { + virtual inline void reinit(const libMesh::Elem& e) { libmesh_assert_msg(_initialized, "Object not initialized."); - context_type c; - c.elem = &e; - c.sys = _sys; - c.physics = _physics; - c.sol = _X; + _context.elem = &e; + _context.sys = _sys; + _context.physics = _physics; + _context.sol = _X; + + _fe_derived->set_compute_dphi_dx(true); + _fe_derived->set_compute_dphi_dx(true); + _fe_derived->set_compute_detJxW(true); + _fe_var_data->set_compute_du_dx(true); _fe_basis->reinit(e); - _fe_derived->reinit(c); - _fe_var_data->init(c); + _fe_derived->reinit(_context); + _fe_var_data->init(_context, {1, 2}); + } + + virtual inline void resize_residual(vector_type& res) const { + + libmesh_assert_msg(_initialized, "Object not initialized"); + res.setZero(_strain_energy->n_dofs()); + } + + virtual inline void resize_jacobian(matrix_type& jac) const { + + libmesh_assert_msg(_initialized, "Object not initialized"); + jac.setZero(_strain_energy->n_dofs(), _strain_energy->n_dofs()); } + + virtual inline void compute(vector_type& res, + matrix_type* jac = nullptr) { + + libmesh_assert_msg(_initialized, "Object not initialized"); + + (*_strain_energy)(_context, res, jac); + } + + virtual inline void compute_sensitivity(const MAST::FunctionBase& f, + vector_type& res) {} + + virtual inline void compute_complex_step_jacobian() {} - virtual void compute(EigenVector& res, - EigenMatrix* jac) {} - virtual void compute_sensitivity(const MAST::FunctionBase& f, - EigenVector& res) {} - virtual void compute_complex_step_jacobian() {} - virtual void compute_auto_diff_jacobian() {} + virtual inline void compute_auto_diff_jacobian() {} protected: @@ -180,6 +209,8 @@ public MAST::ComputationBase { fe_derived_type *_fe_derived; fe_var_type *_fe_var_data; typename MAST::Element2D::LinearContinuumStrainEnergy *_strain_energy; + section_property_type *_section_property; + context_type _context; }; } diff --git a/src/elasticity/strain_energy_compute_kernel.h b/src/elasticity/strain_energy_compute_kernel.h index 888416fb..56f2363c 100644 --- a/src/elasticity/strain_energy_compute_kernel.h +++ b/src/elasticity/strain_energy_compute_kernel.h @@ -168,8 +168,17 @@ public ComputeKernel { _fe_var_data = &fe_data; } + virtual inline uint_type n_dofs() const { + + libmesh_assert_msg(_fe_var_data, "FE data not initialized."); + return 2*_fe_var_data->get_fe_shape_data().n_basis(); + } + virtual inline void operator() (ContextType& c, vector_type& res, matrix_type* jac = nullptr) const { + libmesh_assert_msg(_fe_var_data, "FE data not initialized."); + libmesh_assert_msg(_property, "Section property not initialized"); + const typename fe_var_type::fe_shape_data_type &fe = _fe_var_data->get_fe_shape_data(); diff --git a/src/mesh/fe.h b/src/mesh/fe.h index 19d18a91..82b16c87 100644 --- a/src/mesh/fe.h +++ b/src/mesh/fe.h @@ -22,8 +22,11 @@ class FEBasis: ComputeKernelBase { { } virtual ~FEBasis() {} - virtual inline void set_quadrature(const MAST::Quadrature& q) - { _quadrature = &q;} + virtual inline void set_quadrature(const MAST::Quadrature& q) { + + libmesh_assert_msg(!_quadrature, "Quadrature already initialized."); + _quadrature = &q; + } virtual inline uint_type dim() const = 0; virtual inline uint_type order() const = 0; @@ -121,6 +124,7 @@ class libMeshFE: public FEBasis { libmesh_assert_msg(_fe, "Object not initialized"); + _fe->get_phi(); if (_compute_dphi_dxi) _fe->get_JxW(); diff --git a/src/mesh/fe_geometry_basis_derived_data.h b/src/mesh/fe_geometry_basis_derived_data.h index 105ee4e2..052c7f4e 100644 --- a/src/mesh/fe_geometry_basis_derived_data.h +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -44,10 +44,12 @@ public MAST::FEShapeDataBasen_q_points(); const int_type - d = this->spatial_dim(); + d = c.elem->dim(); const uint_type n_nodes = c.n_nodes(); + this->_spatial_dim = d; + // for this class the number of basis functions should be equal to the number // of nodes libmesh_assert_equal_to(n_nodes, this->_fe_basis->n_basis()); diff --git a/src/mesh/fe_shape_data_base.h b/src/mesh/fe_shape_data_base.h index e12192ce..80aae8f9 100644 --- a/src/mesh/fe_shape_data_base.h +++ b/src/mesh/fe_shape_data_base.h @@ -67,7 +67,7 @@ public MAST::ComputeKernelBase { inline void set_compute_detJxW(bool f) { _compute_JxW = f; - if (f) this->set_compute_Jac(true); + if (f) this->set_compute_detJ(true); } inline void set_compute_dphi_dx(bool f) { @@ -178,7 +178,11 @@ public MAST::ComputeKernelBase { } virtual inline NodalScalarType dphi_dx(uint_type qp, uint_type phi_i, uint_type x_i) const - { return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i);} + { + libmesh_assert_msg(_compute_dphi_dx, "Jacobian inverse computation not requested"); + + return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i); + } protected: diff --git a/src/mesh/fe_var_data.h b/src/mesh/fe_var_data.h index ccbc7044..af0f3387 100644 --- a/src/mesh/fe_var_data.h +++ b/src/mesh/fe_var_data.h @@ -30,7 +30,7 @@ public MAST::ComputeKernelBase { using var_scalar_type = typename MAST::DeducedScalarType::type; using vector_type = typename EigenTraits::vector_type; - using matrix_type = typename EigenTraits::vector_type; + using matrix_type = typename EigenTraits::matrix_type; using fe_shape_data_type = typename MAST::FEShapeDataBase; FEVarData(const std::string& nm): @@ -47,67 +47,72 @@ public MAST::ComputeKernelBase { libmesh_assert_msg(!_coeffs, "Object already initialized"); libmesh_assert_equal_to(_fe->spatial_dim(), c.elem->dim()); + libmesh_assert_less_equal(active_vars.size(), c.sys->n_vars()); + uint_type + d = _fe->spatial_dim(), + n_qp = _fe->n_q_points(), + n_basis = _fe->n_basis(), + n_comp = active_vars.size(); + + _active_vars = active_vars; + // first get the coefficients. If the vector coeffs is provided then // use that, else extract it from the context. if (coeffs) { + _coeffs = coeffs; _coeff_vec.setZero(); } else { std::vector dofs; - c.sys->get_dof_map().dof_indices(c.elem, dofs); - - uint_type - n_dofs = dofs.size(); - - _coeff_vec = vector_type::Zero(n_dofs); - - for (uint_type i=0; iel(dofs[i]); + + _coeff_vec = vector_type::Zero(n_comp*n_basis); + + for (uint_type i=0; iget_dof_map().dof_indices(c.elem, dofs, active_vars[i]); + + for (uint_type j=0; jel(dofs[j]); + } _coeffs = &_coeff_vec; } - // number of avtive variables should be less than or equal to the number of - // variables in the system - libmesh_assert_less_equal(active_vars.size(), c.sys->n_vars()); - - _active_vars = active_vars; - - uint_type - d = _fe->spatial_dim(), - n_qp = _fe->n_q_points(), - n_basis = _fe->n_basis(), - n_comp = active_vars.size(); - libmesh_assert_equal_to(_coeffs->size(), n_comp*n_basis); - _u = vector_type::Zero(n_qp, n_comp); + _u = matrix_type::Zero(n_qp, n_comp); // now, initialize the solution value and derivatives. for (uint_type i=0; iphi(i, j) * (*_coeffs)(j*n_comp+k); + _u(i, j) += _fe->phi(i, j) * (*_coeffs)(_active_vars[j]*n_comp+k); if (_compute_du_dx) { - _du_dx = vector_type::Zero(n_qp, n_comp*d); + _du_dx = matrix_type::Zero(n_qp, n_comp*d); for (uint_type i=0; idphi_dx(i, k, l) * (*_coeffs)(j*n_comp+k); + _du_dx(i, d*j+l) += _fe->dphi_dx(i, k, l) * (*_coeffs)(_active_vars[j]*n_comp+k); } + else + _du_dx.setZero(); } virtual inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - virtual inline void set_fe_shape_data(const fe_shape_data_type& fe) { _fe = &fe;} + virtual inline void set_fe_shape_data(const fe_shape_data_type& fe) { + + libmesh_assert_msg(!_fe, "FE pointer already initialized."); + _fe = &fe; + } virtual inline const fe_shape_data_type& get_fe_shape_data() const { @@ -131,6 +136,7 @@ public MAST::ComputeKernelBase { virtual inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) const { + libmesh_assert_msg(_compute_du_dx, "Object not initialized with du/dx"); libmesh_assert_msg(_coeffs, "Object not initialized"); libmesh_assert_less_equal(comp, _active_vars.size()); diff --git a/src/numerics/fem_operator_matrix.h b/src/numerics/fem_operator_matrix.h index 7bd0ad04..5021756c 100644 --- a/src/numerics/fem_operator_matrix.h +++ b/src/numerics/fem_operator_matrix.h @@ -270,10 +270,11 @@ set_shape_function(unsigned int interpolated_var, ScalarType* vec = _var_shape_functions[discrete_var*_n_interpolated_vars+interpolated_var]; - libmesh_assert_msg(!vec, "Shape function already set"); - - vec = new ScalarType[shape_func.size()]; - _var_shape_functions[discrete_var*_n_interpolated_vars+interpolated_var] = vec; + if (!vec) { + + vec = new ScalarType[shape_func.size()]; + _var_shape_functions[discrete_var*_n_interpolated_vars+interpolated_var] = vec; + } for (uint_type i=0; i<_n_dofs_per_var; i++) vec[i] = shape_func(i); From 11bd7b0a438e024865674d5857afc229c419b1b6 Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Thu, 11 Jun 2020 22:26:37 -0500 Subject: [PATCH 9/9] -- multiple bug-fixes for the implementation of finite element quadrature and derivative computations. -- implementation of complex-step derivative for computation of element Jacobian. --- examples/structural/example_6/example_6.cpp | 90 +++++++++++++- src/base/mast_data_types.h | 1 + src/elasticity/elasticity_elem_operations.h | 67 ++++++++-- src/mesh/fe_geometry_basis_derived_data.h | 2 +- src/mesh/fe_shape_data_base.h | 2 +- src/mesh/fe_var_data.h | 131 ++++++++++++++------ 6 files changed, 241 insertions(+), 52 deletions(-) diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp index 9e017ae3..270899b2 100644 --- a/examples/structural/example_6/example_6.cpp +++ b/examples/structural/example_6/example_6.cpp @@ -1362,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"); @@ -1403,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); @@ -1449,11 +1449,89 @@ int main(int argc, char* argv[]) { else optimizer->optimize(); } + */ - MAST::ElasticityElemOperations> elem_ops; - MAST::ElasticityElemOperations> elem_ops_complex_shape; - MAST::ElasticityElemOperations> elem_ops_complex_sol; - MAST::ElasticityElemOperations> elem_ops_complex_sol2; + TopologyOptimizationSIMP + *o = dynamic_cast*>(top_opt.get()); + const libMesh::Elem* e = + *(o->_mesh->elements_begin()); + e->print_info(); + + EigenMatrix::type jac0; + { + MAST::ElasticityElemOperations> elem_ops; + + elem_ops.init(*o->_sys, *o->_discipline, *o->_sys->solution); + elem_ops.reinit(*e); + + EigenVector::type v; + EigenMatrix::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> elem_ops_complex_shape; + + elem_ops_complex_shape.init(*o->_sys, *o->_discipline, *o->_sys->solution); + elem_ops_complex_shape.reinit(*e); + + EigenVector::type v; + EigenMatrix::type jac; + EigenMatrix::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(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> elem_ops_complex_sol; + + elem_ops_complex_sol.init(*o->_sys, *o->_discipline, *o->_sys->solution); + elem_ops_complex_sol.reinit(*e); + + EigenVector::type v; + EigenMatrix::type jac; + EigenMatrix::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(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> elem_ops_complex_sol2; + + elem_ops_complex_sol2.init(*o->_sys, *o->_discipline, *o->_sys->solution); + elem_ops_complex_sol2.reinit(*e); + + EigenVector::type v; + EigenMatrix::type jac; + EigenMatrix::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(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; diff --git a/src/base/mast_data_types.h b/src/base/mast_data_types.h index ae020c41..793b97b7 100644 --- a/src/base/mast_data_types.h +++ b/src/base/mast_data_types.h @@ -33,6 +33,7 @@ typedef unsigned int uint_type; typedef int int_type; typedef libMesh::Real Real; typedef libMesh::Complex Complex; +#define ComplexStepDelta 1.e-12 typedef Matrix RealVectorX; typedef Matrix RealVector3; diff --git a/src/elasticity/elasticity_elem_operations.h b/src/elasticity/elasticity_elem_operations.h index 18aaa9b0..df390f46 100644 --- a/src/elasticity/elasticity_elem_operations.h +++ b/src/elasticity/elasticity_elem_operations.h @@ -55,6 +55,15 @@ struct ElasticityElemOperationsContext { physics (nullptr), sol (nullptr) {} + void clear() { + + qp = 0; + elem = nullptr; + sys = nullptr; + physics = nullptr; + sol = nullptr; + } + uint_type qp; const libMesh::Elem* elem; const MAST::NonlinearSystem* sys; @@ -116,6 +125,10 @@ public MAST::ComputationBase { _sys = &sys; _physics = &physics; _X = / + + _context.sys = _sys; + _context.physics = _physics; + _context.sol = _X; _quadrature = new quadrature_type("quadrature"); _fe_basis = new fe_basis_type("fe_basis"); @@ -143,6 +156,7 @@ public MAST::ComputationBase { _sys = nullptr; _physics = nullptr; + _context.clear(); delete _quadrature; _quadrature = nullptr; delete _fe_basis; _fe_basis = nullptr; delete _fe_derived; _fe_derived = nullptr; @@ -154,12 +168,9 @@ public MAST::ComputationBase { virtual inline void reinit(const libMesh::Elem& e) { libmesh_assert_msg(_initialized, "Object not initialized."); - + _context.elem = &e; - _context.sys = _sys; - _context.physics = _physics; - _context.sol = _X; - + _fe_derived->set_compute_dphi_dx(true); _fe_derived->set_compute_dphi_dx(true); _fe_derived->set_compute_detJxW(true); @@ -170,13 +181,15 @@ public MAST::ComputationBase { _fe_var_data->init(_context, {1, 2}); } - virtual inline void resize_residual(vector_type& res) const { + template + inline void resize_residual(VecType& res) const { libmesh_assert_msg(_initialized, "Object not initialized"); res.setZero(_strain_energy->n_dofs()); } - virtual inline void resize_jacobian(matrix_type& jac) const { + template + inline void resize_jacobian(MatType& jac) const { libmesh_assert_msg(_initialized, "Object not initialized"); jac.setZero(_strain_energy->n_dofs(), _strain_energy->n_dofs()); @@ -192,13 +205,49 @@ public MAST::ComputationBase { virtual inline void compute_sensitivity(const MAST::FunctionBase& f, vector_type& res) {} + + template + inline void compute_complex_step_jacobian(typename EigenMatrix::type& jac) { + + libmesh_assert_msg(false, "Method only called with Real-valued matrices"); + } - virtual inline void compute_complex_step_jacobian() {} + + template <> + inline void compute_complex_step_jacobian(typename EigenMatrix::type& jac) { + + libmesh_assert_msg(_initialized, "Object not initialized."); + libmesh_assert_msg(_context.elem, "reinit() not called before complex-step computation."); + + vector_type res; + this->resize_residual(res); + this->resize_jacobian(jac); + + for (uint_type i=0; iclear_coeffs_and_vars(); + this->resize_residual(res); + this->_reinit_for_complex_step_jacobian(i); + this->compute(res); + jac.col(i) = res.imag(); + } + + jac /= (ComplexStepDelta); + } virtual inline void compute_auto_diff_jacobian() {} protected: - + + virtual inline void _reinit_for_complex_step_jacobian(uint_type idx) { + + libmesh_assert_msg(_initialized, "Object not initialized."); + libmesh_assert_msg(_context.elem, "reinit() not called before complex-step computation."); + + _fe_var_data->init_for_complex_step_sens(_context, idx, {1, 2}); + } + + bool _initialized; const MAST::NonlinearSystem* _sys; const MAST::PhysicsDisciplineBase* _physics; diff --git a/src/mesh/fe_geometry_basis_derived_data.h b/src/mesh/fe_geometry_basis_derived_data.h index 052c7f4e..550c4bb2 100644 --- a/src/mesh/fe_geometry_basis_derived_data.h +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -139,7 +139,7 @@ public MAST::FEShapeDataBase_compute_dphi_dx) { - this->_dphi_dx = EigenMatrix::type::Zero(nq, d*d*this->n_basis()); + this->_dphi_dx = EigenMatrix::type::Zero(nq, d*this->n_basis()); for (uint_type i=0; i { { libmesh_assert_msg(_compute_dphi_dx, "Jacobian inverse computation not requested"); - return _dphi_dx(qp, x_i*this->spatial_dim()+phi_i); + return _dphi_dx(qp, x_i*this->n_basis()+phi_i); } protected: diff --git a/src/mesh/fe_var_data.h b/src/mesh/fe_var_data.h index af0f3387..47264be1 100644 --- a/src/mesh/fe_var_data.h +++ b/src/mesh/fe_var_data.h @@ -44,6 +44,78 @@ public MAST::ComputeKernelBase { virtual inline void init(const ContextType& c, const std::vector& active_vars, const vector_type* coeffs = nullptr) { + + this->_init_coefficients(c, active_vars, coeffs); + this->_init_variables(c); + } + + virtual inline void init_for_complex_step_sens(const ContextType& c, + const uint_type complex_step_dof, + const std::vector& active_vars) { + + this->_init_coefficients(c, active_vars, nullptr); + + // the provided coefficients shoudl not already have + libmesh_assert_equal_to(_coeffs->imag().norm(), 0.); + + // add complex-step perturbation to the complex-step dof specified + libmesh_assert_less_equal(complex_step_dof, _coeffs->size()); + this->_add_complex_perturbation(_coeff_vec, complex_step_dof); + + this->_init_variables(c); + } + + inline void clear_coeffs_and_vars() { + + _coeffs = nullptr; + _coeff_vec.setZero(); + + _u.setZero(); + _du_dx.setZero(); + } + + virtual inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} + + virtual inline void set_fe_shape_data(const fe_shape_data_type& fe) { + + libmesh_assert_msg(!_fe, "FE pointer already initialized."); + _fe = &fe; + } + + virtual inline const fe_shape_data_type& get_fe_shape_data() const { + + libmesh_assert_msg(_fe, "FE pointer not initialized."); + return *_fe; + } + + virtual inline uint_type n_q_points() const { + + libmesh_assert_msg(_fe, "FE pointer not initialized"); + return _fe->n_q_points(); + } + + virtual inline var_scalar_type u(uint_type qp, uint_type comp) const { + + libmesh_assert_msg(_coeffs, "Object not initialized"); + libmesh_assert_less_equal(comp, _active_vars.size()); + + return _u(qp, comp); + } + + virtual inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) const + { + libmesh_assert_msg(_compute_du_dx, "Object not initialized with du/dx"); + libmesh_assert_msg(_coeffs, "Object not initialized"); + libmesh_assert_less_equal(comp, _active_vars.size()); + + return _du_dx(qp, _fe->spatial_dim()*comp + x_i); + } + +protected: + + inline void _init_coefficients(const ContextType& c, + const std::vector& active_vars, + const vector_type* coeffs = nullptr) { libmesh_assert_msg(!_coeffs, "Object already initialized"); libmesh_assert_equal_to(_fe->spatial_dim(), c.elem->dim()); @@ -82,6 +154,20 @@ public MAST::ComputeKernelBase { } libmesh_assert_equal_to(_coeffs->size(), n_comp*n_basis); + } + + + inline void _init_variables(const ContextType& c) { + + libmesh_assert_msg(_coeffs, "Coefficients not initialized"); + libmesh_assert_equal_to(_fe->spatial_dim(), c.elem->dim()); + libmesh_assert_less_equal(_active_vars.size(), c.sys->n_vars()); + + uint_type + d = _fe->spatial_dim(), + n_qp = _fe->n_q_points(), + n_basis = _fe->n_basis(), + n_comp = _active_vars.size(); _u = matrix_type::Zero(n_qp, n_comp); @@ -89,7 +175,7 @@ public MAST::ComputeKernelBase { for (uint_type i=0; iphi(i, j) * (*_coeffs)(_active_vars[j]*n_comp+k); + _u(i, j) += _fe->phi(i, k) * (*_coeffs)(j*n_basis+k); if (_compute_du_dx) { @@ -99,52 +185,27 @@ public MAST::ComputeKernelBase { for (uint_type j=0; jdphi_dx(i, k, l) * (*_coeffs)(_active_vars[j]*n_comp+k); + _du_dx(i, d*j+l) += _fe->dphi_dx(i, k, l) * (*_coeffs)(j*n_basis+k); } else _du_dx.setZero(); } - virtual inline void set_compute_du_dx(bool f) { _compute_du_dx = f;} - - virtual inline void set_fe_shape_data(const fe_shape_data_type& fe) { - - libmesh_assert_msg(!_fe, "FE pointer already initialized."); - _fe = &fe; - } - - virtual inline const fe_shape_data_type& get_fe_shape_data() const { + template + inline void _add_complex_perturbation(vector_type& v, uint_type i) { - libmesh_assert_msg(_fe, "FE pointer not initialized."); - return *_fe; - } - - virtual inline uint_type n_q_points() const { + libmesh_assert_less_equal(i, v.size()); - libmesh_assert_msg(_fe, "FE pointer not initialized"); - return _fe->n_q_points(); + v(i) += Complex(0, ComplexStepDelta); } - - virtual inline var_scalar_type u(uint_type qp, uint_type comp) const { - - libmesh_assert_msg(_coeffs, "Object not initialized"); - libmesh_assert_less_equal(comp, _active_vars.size()); - return _u(qp, comp); - } - - virtual inline var_scalar_type du_dx(uint_type qp, uint_type comp, uint_type x_i) const - { - libmesh_assert_msg(_compute_du_dx, "Object not initialized with du/dx"); - libmesh_assert_msg(_coeffs, "Object not initialized"); - libmesh_assert_less_equal(comp, _active_vars.size()); - - return _du_dx(qp, _fe->spatial_dim()*comp + x_i); + template <> + inline void _add_complex_perturbation(vector_type& v, uint_type i) { + + libmesh_assert_msg(false, "Complex perturbation cannot be added to real vectors"); } -protected: - bool _compute_du_dx; const fe_shape_data_type *_fe; const vector_type *_coeffs;