meshmover_linear_elasticity.hpp

#ifndef __MESHMOVER_LINEAR_ELASTICITY_H__
#define __MESHMOVER_LINEAR_ELASTICITY_H__

#include <deal.II/lac/trilinos_sparse_matrix.h>

#include "parameters/all_parameters.h"

#include "high_order_grid.h"

namespace PHiLiP {

namespace MeshMover
{
    template <int dim = PHILIP_DIM, typename real = double>
    class LinearElasticity
    {
    using VectorType = dealii::LinearAlgebra::distributed::Vector<real>;
    using DoFHandlerType = dealii::DoFHandler<dim>;
#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D

    using Triangulation = dealii::Triangulation<dim>;
#else

    using Triangulation = dealii::parallel::distributed::Triangulation<dim>;
#endif
      public:
        LinearElasticity(
            const Triangulation &_triangulation,
            const std::shared_ptr<dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>> mapping_fe_field,
            const DoFHandlerType &_dof_handler,
            const dealii::LinearAlgebra::distributed::Vector<int> &_boundary_ids_vector,
            const dealii::LinearAlgebra::distributed::Vector<double> &_boundary_displacements_vector);

        LinearElasticity(
            const HighOrderGrid<dim,real> &high_order_grid,
   const dealii::LinearAlgebra::distributed::Vector<double> &boundary_displacements_vector);

        VectorType get_volume_displacements();

  void evaluate_dXvdXs();

        void
        apply_dXvdXvs(std::vector<dealii::LinearAlgebra::distributed::Vector<double>> &list_of_vectors, dealii::TrilinosWrappers::SparseMatrix &output_matrix);

        void
        apply_dXvdXvs(const dealii::LinearAlgebra::distributed::Vector<double> &input_vector, dealii::LinearAlgebra::distributed::Vector<double> &output_vector);

        void
        apply_dXvdXvs_transpose(
            const dealii::LinearAlgebra::distributed::Vector<double> &input_vector,
            dealii::LinearAlgebra::distributed::Vector<double> &output_vector);

        VectorType displacement_solution;
        dealii::AffineConstraints<double> hanging_node_constraints;

        //dealii::TrilinosWrappers::SparseMatrix dXvdXs;
        std::vector<dealii::LinearAlgebra::distributed::Vector<double>> dXvdXs;

        // /** Sparse matrix containing the dXvdXs sensititivies.
        //  */
        // dealii::TrilinosWrappers::SparseMatrix<double> dXvdXs_matrix;

      private:
        void setup_system();
        void assemble_system();


        void solve_timestep();

        unsigned int solve_linear_problem();

        const Triangulation &triangulation; 
        const std::shared_ptr<dealii::MappingFEField<dim,dim,VectorType,DoFHandlerType>> mapping_fe_field;
        const DoFHandlerType &dof_handler;
        const dealii::QGauss<dim> quadrature_formula;

        dealii::TrilinosWrappers::SparseMatrix system_matrix_unconstrained;
        dealii::LinearAlgebra::distributed::Vector<double> system_rhs_unconstrained;

        dealii::TrilinosWrappers::SparseMatrix system_matrix;
        dealii::LinearAlgebra::distributed::Vector<double> system_rhs;

        dealii::AffineConstraints<double> all_constraints;
        dealii::AffineConstraints<double> dirichlet_boundary_constraints;

        MPI_Comm mpi_communicator; 
        const unsigned int n_mpi_processes; 
        const unsigned int this_mpi_process; 
        dealii::ConditionalOStream pcout; 
        std::vector<dealii::types::global_dof_index> local_dofs_per_process; 
        dealii::IndexSet locally_owned_dofs; 
        dealii::IndexSet locally_relevant_dofs; 

        const dealii::LinearAlgebra::distributed::Vector<int> &boundary_ids_vector;
        const dealii::LinearAlgebra::distributed::Vector<double> &boundary_displacements_vector;

        dealii::LinearAlgebra::distributed::Vector<double> tensor_to_vector(const std::vector<dealii::Tensor<1,dim,real>> &boundary_displacements_tensors) const;

    };
} // namespace MeshMover

} // namespace PHiLiP

#endif