Helper function for converting std::vector into `Tensor<VectorizedA…

…rray> ` (#383)

* removing some datatypes

* adding function to convert std::vector to 1st rank tensors of vectorized arrays

include/core/matrixFreePDE.h

@@ -1,4 +1,3 @@
-// base class for matrix Free implementation of PDE's
 #ifndef MATRIXFREEPDE_H
 #define MATRIXFREEPDE_H
 
@@ -36,20 +35,10 @@
 #include <grains/SimplifiedGrainRepresentation.h>
 #include <nucleation/nucleus.h>
 #include <utilities/computeStress.h>
+#include <utilities/utilities.h>
 
 using namespace dealii;
 
-// define data types
-#ifndef scalarType
-using scalarType = VectorizedArray<double>;
-#endif
-#ifndef vectorType
-using vectorType = LinearAlgebra::distributed::Vector<double>;
-#endif
-
-// macro for constants
-#define constV(a) make_vectorized_array(a)
-
 /**
  * \brief This is the abstract base class for the matrix free implementation of parabolic
  * and elliptic BVP's, and supports MPI, threads and vectorization (Hybrid Parallel). This
@@ -118,7 +107,8 @@ class MatrixFreePDE : public Subscriptor
    * \param src The source vector.
    */
   void
-  vmult(vectorType &dst, const vectorType &src) const;
+  vmult(dealii::LinearAlgebra::distributed::Vector<double>       &dst,
+        const dealii::LinearAlgebra::distributed::Vector<double> &src) const;
 
   /**
    * \brief Vector of all the physical fields in the problem. Fields are identified by
@@ -256,13 +246,16 @@ class MatrixFreePDE : public Subscriptor
   std::vector<IndexSet *> locally_relevant_dofsSet_nonconst;
   /*Vector all the solution vectors in the problem. In a multi-field problem, each primal
    * field has a solution vector associated with it.*/
-  std::vector<vectorType *> solutionSet;
+  std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> solutionSet;
   /*Vector all the residual (RHS) vectors in the problem. In a multi-field problem, each
    * primal field has a residual vector associated with it.*/
-  std::vector<vectorType *> residualSet;
+  std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> residualSet;
   /*Vector of parallel solution transfer objects. This is used only when adaptive meshing
    * is enabled.*/
-  std::vector<parallel::distributed::SolutionTransfer<dim, vectorType> *> soltransSet;
+  std::vector<parallel::distributed::SolutionTransfer<
+    dim,
+    dealii::LinearAlgebra::distributed::Vector<double>> *>
+    soltransSet;
 
   // matrix free objects
   /*Object of class MatrixFree<dim>. This is primarily responsible for all the
@@ -273,14 +266,14 @@ class MatrixFreePDE : public Subscriptor
   /*Vector to store the inverse of the mass matrix diagonal for scalar fields.
    * Due to the choice of spectral elements with Guass-Lobatto quadrature, the
    * mass matrix is diagonal.*/
-  vectorType invMscalar;
+  dealii::LinearAlgebra::distributed::Vector<double> invMscalar;
   /*Vector to store the inverse of the mass matrix diagonal for vector fields.
    * Due to the choice of spectral elements with Guass-Lobatto quadrature, the
    * mass matrix is diagonal.*/
-  vectorType invMvector;
+  dealii::LinearAlgebra::distributed::Vector<double> invMvector;
   /*Vector to store the solution increment. This is a temporary vector used
    * during implicit solves of the Elliptic fields.*/
-  vectorType dU_vector, dU_scalar;
+  dealii::LinearAlgebra::distributed::Vector<double> dU_vector, dU_scalar;
 
   // matrix free methods
   /*Current field index*/
@@ -321,42 +314,44 @@ class MatrixFreePDE : public Subscriptor
   // virtual methods to be implemented in the derived class
   /*Method to calculate LHS(implicit solve)*/
   void
-  getLHS(const MatrixFree<dim, double>               &data,
-         vectorType                                  &dst,
-         const vectorType                            &src,
-         const std::pair<unsigned int, unsigned int> &cell_range) const;
+  getLHS(const MatrixFree<dim, double>                            &data,
+         dealii::LinearAlgebra::distributed::Vector<double>       &dst,
+         const dealii::LinearAlgebra::distributed::Vector<double> &src,
+         const std::pair<unsigned int, unsigned int>              &cell_range) const;
 
   bool generatingInitialGuess;
   void
-  getLaplaceLHS(const MatrixFree<dim, double>               &data,
-                vectorType                                  &dst,
-                const vectorType                            &src,
+  getLaplaceLHS(const MatrixFree<dim, double>                            &data,
+                dealii::LinearAlgebra::distributed::Vector<double>       &dst,
+                const dealii::LinearAlgebra::distributed::Vector<double> &src,
                 const std::pair<unsigned int, unsigned int> &cell_range) const;
 
   void
   setNonlinearEqInitialGuess();
   void
   computeLaplaceRHS(unsigned int fieldIndex);
   void
-  getLaplaceRHS(const MatrixFree<dim, double>               &data,
-                vectorType                                  &dst,
-                const vectorType                            &src,
+  getLaplaceRHS(const MatrixFree<dim, double>                            &data,
+                dealii::LinearAlgebra::distributed::Vector<double>       &dst,
+                const dealii::LinearAlgebra::distributed::Vector<double> &src,
                 const std::pair<unsigned int, unsigned int> &cell_range) const;
 
   /*Method to calculate RHS (implicit/explicit). This is an abstract method, so
    * every model which inherits MatrixFreePDE<dim> has to implement this
    * method.*/
   void
-  getExplicitRHS(const MatrixFree<dim, double>               &data,
-                 std::vector<vectorType *>                   &dst,
-                 const std::vector<vectorType *>             &src,
-                 const std::pair<unsigned int, unsigned int> &cell_range) const;
+  getExplicitRHS(
+    const MatrixFree<dim, double>                                           &data,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *>       &dst,
+    const std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &src,
+    const std::pair<unsigned int, unsigned int> &cell_range) const;
 
   void
-  getNonexplicitRHS(const MatrixFree<dim, double>               &data,
-                    std::vector<vectorType *>                   &dst,
-                    const std::vector<vectorType *>             &src,
-                    const std::pair<unsigned int, unsigned int> &cell_range) const;
+  getNonexplicitRHS(
+    const MatrixFree<dim, double>                                           &data,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *>       &dst,
+    const std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &src,
+    const std::pair<unsigned int, unsigned int> &cell_range) const;
 
   virtual void
   explicitEquationRHS(
@@ -387,13 +382,15 @@ class MatrixFreePDE : public Subscriptor
     [[maybe_unused]] const Point<dim, VectorizedArray<double>> q_point_loc,
     [[maybe_unused]] const VectorizedArray<double>             element_volume) const {};
   void
-  computePostProcessedFields(std::vector<vectorType *> &postProcessedSet);
+  computePostProcessedFields(
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &postProcessedSet);
 
   void
-  getPostProcessedFields(const MatrixFree<dim, double>               &data,
-                         std::vector<vectorType *>                   &dst,
-                         const std::vector<vectorType *>             &src,
-                         const std::pair<unsigned int, unsigned int> &cell_range);
+  getPostProcessedFields(
+    const MatrixFree<dim, double>                                           &data,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *>       &dst,
+    const std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &src,
+    const std::pair<unsigned int, unsigned int>                             &cell_range);
 
   // methods to apply dirichlet BC's
   /*Map of degrees of freedom to the corresponding Dirichlet boundary
@@ -506,9 +503,10 @@ class MatrixFreePDE : public Subscriptor
 
   /*Method to compute the integral of a field.*/
   void
-  computeIntegral(double                   &integratedField,
-                  int                       index,
-                  std::vector<vectorType *> variableSet);
+  computeIntegral(
+    double                                                           &integratedField,
+    int                                                               index,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> variableSet);
 
   // variables for time dependent problems
   /*Flag used to see if invM, time stepping in run(), etc are necessary*/

include/core/refinement/AdaptiveRefinement.h

@@ -22,18 +22,18 @@ template <int dim, int degree>
 class AdaptiveRefinement
 {
 public:
-  using vectorType = dealii::LinearAlgebra::distributed::Vector<double>;
-
   /**
    * Default constructor.
    */
   AdaptiveRefinement(
-    const userInputParameters<dim>            &_userInputs,
-    parallel::distributed::Triangulation<dim> &_triangulation,
-    std::vector<Field<dim>>                   &_fields,
-    std::vector<vectorType *>                 &_solutionSet,
-    std::vector<parallel::distributed::SolutionTransfer<dim, vectorType> *> &_soltransSet,
-    std::vector<FESystem<dim> *>                                            &_FESet,
+    const userInputParameters<dim>                                    &_userInputs,
+    parallel::distributed::Triangulation<dim>                         &_triangulation,
+    std::vector<Field<dim>>                                           &_fields,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &_solutionSet,
+    std::vector<parallel::distributed::SolutionTransfer<
+      dim,
+      dealii::LinearAlgebra::distributed::Vector<double>> *>          &_soltransSet,
+    std::vector<FESystem<dim> *>                                      &_FESet,
     std::vector<DoFHandler<dim> *>                 &_dofHandlersSet_nonconst,
     std::vector<const AffineConstraints<double> *> &_constraintsDirichletSet,
     std::vector<const AffineConstraints<double> *> &_constraintsOtherSet);
@@ -65,9 +65,11 @@ class AdaptiveRefinement
 
   std::vector<Field<dim>> &fields;
 
-  std::vector<vectorType *> &solutionSet;
+  std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &solutionSet;
 
-  std::vector<parallel::distributed::SolutionTransfer<dim, vectorType> *> &soltransSet;
+  std::vector<parallel::distributed::SolutionTransfer<
+    dim,
+    dealii::LinearAlgebra::distributed::Vector<double>> *> &soltransSet;
 
   std::vector<FESystem<dim> *> &FESet;
 

include/core/typeDefs.h

@@ -8,36 +8,6 @@
 // #ifndef INCLUDE_TYPEDEFS_H_
 // #define INCLUDE_TYPEDEFS_H_
 
-// #include <deal.II/base/quadrature.h>
-// #include <deal.II/base/timer.h>
-// #include <deal.II/lac/vector.h>
-// #include <deal.II/lac/affine_constraints.h>
-// #include <deal.II/fe/fe_system.h>
-// #include <deal.II/fe/fe_q.h>
-// #include <deal.II/fe/fe_values.h>
-// #include <deal.II/grid/tria.h>
-// #include <deal.II/grid/tria_accessor.h>
-// #include <deal.II/grid/tria_iterator.h>
-// #include <deal.II/grid/grid_tools.h>
-// #include <deal.II/dofs/dof_tools.h>
-// #include <deal.II/dofs/dof_handler.h>
-//  #include <deal.II/numerics/vector_tools.h>
-//  #include <deal.II/lac/la_parallel_vector.h>
-//  #include <deal.II/matrix_free/matrix_free.h>
-//  #include <deal.II/matrix_free/fe_evaluation.h>
-//  #include <deal.II/base/config.h>
-//  #include <deal.II/base/exceptions.h>
-//  #include <deal.II/distributed/tria.h>
-// #include <deal.II/distributed/solution_transfer.h>
-// #include <deal.II/grid/manifold_lib.h>
-
-// define data types
-#ifndef scalarType
-using scalarType = dealii::VectorizedArray<double>;
-#endif
-#ifndef vectorType
-using vectorType = dealii::LinearAlgebra::distributed::Vector<double>;
-#endif
 // define FE system types
 #ifndef typeScalar
 using typeScalar = dealii::FEEvaluation<dim, degree, degree + 1, 1, double>;

include/core/variableContainer.h

@@ -89,11 +89,14 @@ class variableContainer
 
   // Initialize, read DOFs, and set evaulation flags for each variable
   void
-  reinit_and_eval(const std::vector<vectorType *> &src, unsigned int cell);
+  reinit_and_eval(
+    const std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &src,
+    unsigned int                                                             cell);
   void
-  reinit_and_eval_change_in_solution(const vectorType &src,
-                                     unsigned int      cell,
-                                     unsigned int      var_being_solved);
+  reinit_and_eval_change_in_solution(
+    const dealii::LinearAlgebra::distributed::Vector<double> &src,
+    unsigned int                                              cell,
+    unsigned int                                              var_being_solved);
 
   // Only initialize the FEEvaluation object for each variable (used for
   // post-processing)
@@ -102,10 +105,12 @@ class variableContainer
 
   // Integrate the residuals and distribute from local to global
   void
-  integrate_and_distribute(std::vector<vectorType *> &dst);
+  integrate_and_distribute(
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &dst);
   void
-  integrate_and_distribute_change_in_solution_LHS(vectorType        &dst,
-                                                  const unsigned int var_being_solved);
+  integrate_and_distribute_change_in_solution_LHS(
+    dealii::LinearAlgebra::distributed::Vector<double> &dst,
+    const unsigned int                                  var_being_solved);
 
   // The quadrature point index, a method to get the number of quadrature points
   // per cell, and a method to get the xyz coordinates for the quadrature point

include/grains/FloodFiller.h

@@ -11,10 +11,6 @@
 #include <deal.II/lac/la_parallel_vector.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
 
-#ifndef vectorType
-using vectorType = dealii::LinearAlgebra::distributed::Vector<double>;
-#endif
-
 /**
  * This class holds information for a grain, including its index and the list of
  * vertices in that grain.
@@ -118,30 +114,30 @@ class FloodFiller
    * mesh/field and outputs a vector of GrainSet objects.
    */
   void
-  calcGrainSets(dealii::FESystem<dim>      &finite_element,
-                dealii::DoFHandler<dim>    &dof_handler,
-                vectorType                 *solution_field,
-                double                      threshold_lower,
-                double                      threshold_upper,
-                int                         min_id,
-                unsigned int                order_parameter_index,
-                std::vector<GrainSet<dim>> &grain_sets);
+  calcGrainSets(dealii::FESystem<dim>                              &finite_element,
+                dealii::DoFHandler<dim>                            &dof_handler,
+                dealii::LinearAlgebra::distributed::Vector<double> *solution_field,
+                double                                              threshold_lower,
+                double                                              threshold_upper,
+                int                                                 min_id,
+                unsigned int                                        order_parameter_index,
+                std::vector<GrainSet<dim>>                         &grain_sets);
 
 protected:
   /**
    * The actual recursive flood fill method.
    */
   template <typename T>
   void
-  recursiveFloodFill(T                           cell,
-                     T                           cell_end,
-                     vectorType                 *solution_field,
-                     double                      threshold_lower,
-                     double                      threshold_upper,
-                     int                         min_id,
-                     unsigned int               &grain_index,
-                     std::vector<GrainSet<dim>> &grain_sets,
-                     bool                       &grain_assigned);
+  recursiveFloodFill(T                                                   cell,
+                     T                                                   cell_end,
+                     dealii::LinearAlgebra::distributed::Vector<double> *solution_field,
+                     double                                              threshold_lower,
+                     double                                              threshold_upper,
+                     int                                                 min_id,
+                     unsigned int                                       &grain_index,
+                     std::vector<GrainSet<dim>>                         &grain_sets,
+                     bool                                               &grain_assigned);
 
   /**
    * The method to merge the grain sets from all the processors.

include/grains/OrderParameterRemapper.h

@@ -21,11 +21,12 @@ class OrderParameterRemapper
    * objects.
    */
   void
-  remap(std::vector<SimplifiedGrainRepresentation<dim>> &grain_representations,
-        std::vector<vectorType *>                       &solution_fields,
-        dealii::DoFHandler<dim>                         &dof_handler,
-        unsigned int                                     dofs_per_cell,
-        double                                           buffer);
+  remap(
+    std::vector<SimplifiedGrainRepresentation<dim>> &grain_representations,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &solution_fields,
+    dealii::DoFHandler<dim>                                           &dof_handler,
+    unsigned int                                                       dofs_per_cell,
+    double                                                             buffer);
 
   /**
    * This method does the core work of the class to reassign grains across
@@ -34,12 +35,12 @@ class OrderParameterRemapper
    */
   void
   remap_from_index_field(
-    std::vector<SimplifiedGrainRepresentation<dim>> &grain_representations,
-    const vectorType                                *grain_index_field,
-    std::vector<vectorType *>                       &solution_fields,
-    dealii::DoFHandler<dim>                         &dof_handler,
-    unsigned int                                     dofs_per_cell,
-    double                                           buffer);
+    std::vector<SimplifiedGrainRepresentation<dim>>          &grain_representations,
+    const dealii::LinearAlgebra::distributed::Vector<double> *grain_index_field,
+    std::vector<dealii::LinearAlgebra::distributed::Vector<double> *> &solution_fields,
+    dealii::DoFHandler<dim>                                           &dof_handler,
+    unsigned int                                                       dofs_per_cell,
+    double                                                             buffer);
 
 protected:
 };