diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp index 3d885464..270899b2 100644 --- a/examples/structural/example_6/example_6.cpp +++ b/examples/structural/example_6/example_6.cpp @@ -51,6 +51,8 @@ #include "examples/structural/base/inplane_2d_model.h" #include "examples/structural/base/truss_2d_model.h" #include "examples/structural/base/eyebar_2d_model.h" +#include "elasticity/elasticity_elem_operations.h" +#include "base/view.h" // libMesh includes @@ -1360,7 +1362,7 @@ _optim_con(int* mode, int main(int argc, char* argv[]) { libMesh::LibMeshInit init(argc, argv); - + MAST::Examples::GetPotWrapper input(argc, argv, "input"); @@ -1401,7 +1403,7 @@ int main(int argc, char* argv[]) { std::string s = input("optimizer", "optimizer to use in the example", "gcmma"); - +/* if (s == "gcmma") { optimizer.reset(new MAST::GCMMAOptimizationInterface); @@ -1447,7 +1449,90 @@ int main(int argc, char* argv[]) { else optimizer->optimize(); } + */ + TopologyOptimizationSIMP + *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/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/computation_base.h b/src/base/computation_base.h new file mode 100644 index 00000000..7bf2f096 --- /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: + + using basis_scalar_type = typename Traits::basis_scalar_type; + using nodal_scalar_type = typename Traits::nodal_scalar_type; + using sol_scalar_type = typename Traits::sol_scalar_type; + + ComputationBase() { } + virtual ~ComputationBase() { } + void add_compute_kernel(MAST::ComputeKernelBase& 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 new file mode 100644 index 00000000..f548f430 --- /dev/null +++ b/src/base/compute_kernel.h @@ -0,0 +1,161 @@ +// +// compute_kernel.h +// compute_kernels +// +// Created by Manav Bhatia on 5/1/20. +// Copyright © 2020 Manav Bhatia. All rights reserved. +// + +#ifndef _mast_compute_kernel_h_ +#define _mast_compute_kernel_h_ + +// C++ includes +#include +#include + +// MAST includes +#include "base/compute_kernel_base.h" +#include "base/view.h" + + +namespace MAST { + + +template +struct IsScalarType { + using type = std::false_type; +}; + + +template <> +struct IsScalarType { + using type = std::true_type; +}; + + +template <> +struct IsScalarType> { + using type = std::true_type; +}; + + + +template ::type> +struct KernelReturnType { + +}; + + +template +struct KernelReturnType { + + using type = KernelValueType; +}; + + +template +struct KernelReturnType { + + using type = KernelValueType; +}; + + + +template +struct KernelReturnType { + + using type = Eigen::Map; +}; + + +template +struct KernelReturnType { + + using type = KernelValueType&; +}; + + +template +KernelViewType build_kernel_view(KernelType& k, ContextType& c) { } + + + + +// Forward decleration +template class ComputeKernelDerivative; + + +template +class ComputeKernel: public MAST::ComputeKernelBase { + +public: + + using derivative_kernel_type = MAST::ComputeKernelDerivative; + + ComputeKernel(const std::string& nm, const bool executable): + MAST::ComputeKernelBase(nm, executable), + _derivative_kernel (nullptr) + {} + + virtual ~ComputeKernel() {} + + virtual inline void set_derivative_kernel(derivative_kernel_type& d) + { + libmesh_assert_msg(!_derivative_kernel, "Derivative kernel already set."); + _derivative_kernel = &d; + } + + virtual inline const derivative_kernel_type& get_derivative_kernel() const + { + libmesh_assert_msg(_derivative_kernel, "Derivative kernel not set"); + return *_derivative_kernel; + } + +protected: + + derivative_kernel_type *_derivative_kernel; +}; + + + +template +class ComputeKernelDerivative: public MAST::ComputeKernelBase { + +public: + + using primal_kernel_type = MAST::ComputeKernel; + + ComputeKernelDerivative (const std::string& nm, const bool executable): + MAST::ComputeKernelBase (nm), + _primal_kernel (nullptr), + _f (nullptr) + {} + + virtual ~ComputeKernelDerivative() {} + + virtual inline void set_primal_kernel(primal_kernel_type& k) + { + libmesh_assert_msg(!_primal_kernel, "Primal kernel already set."); + _primal_kernel = &k; + } + + virtual inline primal_kernel_type& get_primal_kernel() + { + libmesh_assert_msg(_primal_kernel, "Primal kernel not set."); + return *_primal_kernel; + } + + virtual inline void set_derivative_paramter(const MAST::FunctionBase& f); + +protected: + + primal_kernel_type *_primal_kernel; + const MAST::FunctionBase *_f; +}; + + + +} + +#endif /* __mast_compute_kernel_h__ */ + diff --git a/src/base/compute_kernel_base.h b/src/base/compute_kernel_base.h new file mode 100644 index 00000000..6e2bca20 --- /dev/null +++ b/src/base/compute_kernel_base.h @@ -0,0 +1,59 @@ + +#ifndef _mast_compute_kernel_base_h_ +#define _mast_compute_kernel_base_h_ + +// C++ includes +#include +#include + +// MAST includes +#include "base/mast_data_types.h" + + +namespace MAST { + +// Forward declerations +class FunctionBase; + +template +class ComputeKernelBase { + +public: + + ComputeKernelBase(const std::string& nm, + const bool executable): + _nm (nm), + _executable (executable) {} + + virtual ~ComputeKernelBase() {} + + virtual inline bool is_executable() const { return _executable;} + virtual inline bool depends_on(const MAST::ComputeKernelBase& d) const + { return _dependency.count(&d);} + virtual inline const std::set*>& get_dependencies() const + { return _dependency;} + +protected: + + virtual inline void _add_dependency(const MAST::ComputeKernelBase& d) { _dependency.insert(&d);} + + const std::string _nm; + const bool _executable; + 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..793b97b7 100644 --- a/src/base/mast_data_types.h +++ b/src/base/mast_data_types.h @@ -29,8 +29,11 @@ #include "Eigen/Dense" using namespace Eigen; +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; @@ -56,4 +59,59 @@ 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 +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..e2b5e051 --- /dev/null +++ b/src/base/view.h @@ -0,0 +1,240 @@ + +#ifndef _mast_view_h_ +#define _mast_view_h_ + +// C++ includes +#include + +// 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]]; + } + +}; + + + + + +class ViewBase { + +public: + ViewBase(const std::string& nm, const std::vector& sizes): + _name (nm), + _dim (sizes.size()), + _size (sizes) { + + // 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() {} + + inline const std::string& name() const { return _name;} + + 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); + + uint_type n = 0; + for (uint_type i=0; i<_dim; i++) { + libmesh_assert_less(indices[i], _size[i]); + n += indices[i] * _slice_size[i+1]; + } + + libmesh_assert_less(n, _slice_size[0]); + + return n; + } + + inline virtual void init() = 0; + +protected: + + const std::string _name; + const uint_type _dim; + const std::vector _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 _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);} + + 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: + + 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; +}; + +} + +#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..df390f46 --- /dev/null +++ b/src/elasticity/elasticity_elem_operations.h @@ -0,0 +1,266 @@ + +#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 +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 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::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; +}; + + +struct ElasticityElemOperationsContext { + + ElasticityElemOperationsContext(): + qp (0), + elem (nullptr), + sys (nullptr), + 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; + 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(); + } +}; + + +template +class ElasticityElemOperations: +public MAST::ComputationBase { + +public: + + using context_type = typename Traits::context_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; + 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), + _section_property (nullptr) + { } + + virtual ~ElasticityElemOperations() { this->clear();} + + virtual inline void init(const MAST::NonlinearSystem& sys, + const MAST::PhysicsDisciplineBase& physics, + const libMesh::NumericVector& sol) { + + libmesh_assert(!_initialized); + + + // setup the volume compute kernels + _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"); + _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); + _section_property = new section_property_type; + _strain_energy->set_section_property(*_section_property); + + _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 inline void clear() { + + if (_initialized) { + + _sys = nullptr; + _physics = nullptr; + _context.clear(); + 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 inline void reinit(const libMesh::Elem& e) { + + libmesh_assert_msg(_initialized, "Object not initialized."); + + _context.elem = &e; + + _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(_context); + _fe_var_data->init(_context, {1, 2}); + } + + template + inline void resize_residual(VecType& res) const { + + libmesh_assert_msg(_initialized, "Object not initialized"); + res.setZero(_strain_energy->n_dofs()); + } + + 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()); + } + + 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) {} + + template + inline void compute_complex_step_jacobian(typename EigenMatrix::type& jac) { + + libmesh_assert_msg(false, "Method only called with Real-valued matrices"); + } + + + 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; + 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::Element2D::LinearContinuumStrainEnergy *_strain_energy; + section_property_type *_section_property; + context_type _context; +}; +} + +#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..56f2363c --- /dev/null +++ b/src/elasticity/strain_energy_compute_kernel.h @@ -0,0 +1,228 @@ + +#ifndef _mast_strain_energy_compute_kernel_h_ +#define _mast_strain_energy_compute_kernel_h_ + +// MAST includes +#include "mesh/fe_compute_kernel.h" +#include "numerics/fem_operator_matrix.h" + + +namespace MAST { + +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); +} + + +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 ~LinearContinuumStrainEnergy() { } + + virtual inline void + 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) + { + libmesh_assert_msg(!_fe_var_data, "FE data already initialized."); + _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(); + + 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 section_property_type *_property; + const fe_var_type *_fe_var_data; +}; +} +} + +#endif // _mast_strain_energy_compute_kernel_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/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..ccb6d536 --- /dev/null +++ b/src/materials/isotropic_hookean_material.h @@ -0,0 +1,77 @@ + + +#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() { } + + 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 = ν + } + + + void inline operator(const ContextType& c, value_type& m) { + + libmesh_assert_msg(_E && _nu, "Property kernels should be provided before execute()."); + + const scalar_type + E = _E->value(c), + nu = _nu->value(c), + G = E/2./(1.+nu); + + 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); + } + + inline void execute_derivative(const ContextType& c, value_type& m) { + + libmesh_assert_msg(_E && _nu, "Property kernels should be provided before execute()."); + + 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); + } + +protected: + + const MAST::ComputeKernel *_E; + const MAST::ComputeKernel *_nu; +}; + +} + +#endif // _mast_isotropic_hookean_material_h_ + diff --git a/src/mesh/fe.h b/src/mesh/fe.h new file mode 100644 index 00000000..82b16c87 --- /dev/null +++ b/src/mesh/fe.h @@ -0,0 +1,191 @@ + +#ifndef _mast_febasis_h_ +#define _mast_febasis_h_ + +// MAST includes +#include "base/compute_kernel.h" +#include "mesh/quadrature.h" + +namespace MAST { + + +template +class FEBasis: ComputeKernelBase { + +public: + + using basis_scalar_type = ScalarType; + + 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) { + + libmesh_assert_msg(!_quadrature, "Quadrature already initialized."); + _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 + { + 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; + +protected: + + const MAST::Quadrature *_quadrature; +}; + + +/*! serves as a wrapper around libMesh FEBase object to provide */ +template +class libMeshFE: public FEBasis { + +public: + + using scalar_type = ScalarType; + using basis_scalar_type = typename MAST::FEBasis::basis_scalar_type; + + libMeshFE(const std::string& nm): + MAST::FEBasis(nm, false), + _compute_dphi_dxi (false), + _own_pointer (false), + _fe (nullptr) + { } + + virtual ~libMeshFE() { + + if (_own_pointer) + delete _fe; + } + + inline void set_compute_dphi_dxi(bool f) { _compute_dphi_dxi = f;} + + virtual inline void init(libMesh::FEBase& fe) { + + libmesh_assert_msg(!_fe, "Object already initialized"); + + _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"); + + _fe->get_phi(); + 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()}; + else + _dphi_dxi.clear(); + } + + + 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 { + + libmesh_assert_msg(_fe, "FE not initialized."); + + return _fe->get_order(); + } + + 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 + { + 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 + { + libmesh_assert_msg(_compute_dphi_dxi, "FE not initialized with derivatives."); + + return (*_dphi_dxi[xi_i])[phi_i][qp]; + } + +protected: + + bool _compute_dphi_dxi; + bool _own_pointer; + 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..7f394015 --- /dev/null +++ b/src/mesh/fe_basis_derived_data.h @@ -0,0 +1,96 @@ + +#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: + + FEBasisDerivedData(const std::string& nm): + MAST::FEShapeDataBase(nm, false), + _fe_geom(nullptr) + {} + + virtual ~FEBasisDerivedData() {} + +}; + + + +/*! 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) + { + 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 + { + 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 + { + libmesh_assert_msg(_fe_geom, "FE Geometry Basis not initialized."); + + return _fe_geom->detJ(qp); + } + + virtual inline NodalScalarType detJxW(uint_type qp) const override + { + 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 + { + 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 + { + 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; +}; +} + +#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..1bb5c2bf --- /dev/null +++ b/src/mesh/fe_compute_kernel.h @@ -0,0 +1,26 @@ + +#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 { + +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, false) { } + 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..550c4bb2 --- /dev/null +++ b/src/mesh/fe_geometry_basis_derived_data.h @@ -0,0 +1,168 @@ + +#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() {} + +protected: + +}; + + + +/*! 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 void reinit(const ContextType& c) override { + + uint_type + nq = this->n_q_points(); + const int_type + 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()); + + 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*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 { } + + +protected: + +}; + +} + +#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..6ef0386d --- /dev/null +++ b/src/mesh/fe_shape_data_base.h @@ -0,0 +1,210 @@ + +#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, 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, false), + _compute_xyz (false), + _compute_Jac (false), + _compute_Jac_inv (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; + 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_detJ(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."); + + _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 + { + 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 + { + + 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 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 + { + libmesh_assert_msg(_compute_xyz, "Nodal and QPoint locations not requested"); + return _xyz(qp, x_i); + } + + virtual inline NodalScalarType detJ(uint_type qp) const + { + libmesh_assert_msg(_compute_detJ, "Jacobian computation not requested"); + return _detJ(qp); + } + + virtual inline NodalScalarType detJxW(uint_type qp) const + { + 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 + { + 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 + { + 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 + { + 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 + { + libmesh_assert_msg(_compute_dphi_dx, "Jacobian inverse computation not requested"); + + return _dphi_dx(qp, x_i*this->n_basis()+phi_i); + } + +protected: + + bool _compute_xyz; + bool _compute_Jac; + bool _compute_Jac_inv; + bool _compute_detJ; + bool _compute_JxW; + bool _compute_dphi_dx; + bool _compute_normal; + uint_type _spatial_dim; + + MAST::FEBasis *_fe_basis; + + typename EigenMatrix::type _node_coord; + 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; +}; + +} +#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..47264be1 --- /dev/null +++ b/src/mesh/fe_var_data.h @@ -0,0 +1,220 @@ + +#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: + + FEVarData(const std::string& nm): MAST::ComputeKernelBase(nm, false) {} + virtual ~FEVarData() {} + +}; + + + + +template +class FEVarData: +public MAST::ComputeKernelBase { + +public: + + using var_scalar_type = typename MAST::DeducedScalarType::type; + using vector_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): + MAST::ComputeKernelBase(nm, false), + _compute_du_dx (nullptr), + _fe (nullptr), + _coeffs (nullptr) + {} + virtual ~FEVarData() {} + + 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()); + 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; + + _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; + } + + 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); + + // now, initialize the solution value and derivatives. + for (uint_type i=0; iphi(i, k) * (*_coeffs)(j*n_basis+k); + + if (_compute_du_dx) { + + _du_dx = matrix_type::Zero(n_qp, n_comp*d); + + for (uint_type i=0; idphi_dx(i, k, l) * (*_coeffs)(j*n_basis+k); + + } + else + _du_dx.setZero(); + } + + template + inline void _add_complex_perturbation(vector_type& v, uint_type i) { + + libmesh_assert_less_equal(i, v.size()); + + v(i) += Complex(0, ComplexStepDelta); + } + + template <> + inline void _add_complex_perturbation(vector_type& v, uint_type i) { + + libmesh_assert_msg(false, "Complex perturbation cannot be added to real vectors"); + } + + 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; +}; + +} + +#endif // _mast_fe_var_data_h_ diff --git a/src/mesh/quadrature.h b/src/mesh/quadrature.h new file mode 100644 index 00000000..0ca62eb2 --- /dev/null +++ b/src/mesh/quadrature.h @@ -0,0 +1,165 @@ + + +#ifndef _mast_quadrature_h_ +#define _mast_quadrature_h_ + +// MAST includes +#include "base/compute_kernel.h" + +// libMesh includes +#include "libmesh/quadrature.h" +#include "libmesh/enum_quadrature_type.h" + + +namespace MAST { + +template +class Quadrature: public MAST::ComputeKernelBase { + +public: + + using scalar_type = ScalarType; + + 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; + 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 Eigen::Matrix& points, + const typename Eigen::Matrix& 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; + Eigen::Matrix _points; + Eigen::Matrix _weights; +}; + + + +/*! serves as a wrapper around libMesh */ +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, false), + _own_pointer (false), + _q (nullptr) + { } + + virtual ~libMeshQuadrature() { + + if (_q && _own_pointer) + delete _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"); + 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: + + bool _own_pointer; + libMesh::QBase* _q; +}; +} + +#endif //_mast_quadrature_h_ 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..5021756c 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,54 @@ 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; + if (!vec) { + + 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 +327,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 +350,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 +377,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 +406,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 +434,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 +469,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 +497,12 @@ left_multiply_transpose(T& r, const T& m) const { for (unsigned int l=0; l