Merge pull request #213 from asartori86/constraints

allow for different constraint matrices for different matrices

Author: luca-heltai

include/pidomus.h

@@ -321,8 +321,8 @@ class piDoMUS : public ParameterAcceptor
   shared_ptr<FiniteElement<dim, spacedim> >       fe;
   shared_ptr<DoFHandler<dim, spacedim> >          dof_handler;
 
-  ConstraintMatrix                          constraints;
-  ConstraintMatrix                          constraints_dot;
+  std::vector<shared_ptr<ConstraintMatrix> >      train_constraints;
+  ConstraintMatrix                                constraints_dot;
 
   LinearOperator<typename LAC::VectorType> jacobian_preconditioner_op;
   LinearOperator<typename LAC::VectorType> jacobian_preconditioner_op_finer;

include/pidomus_signals.h

@@ -45,10 +45,15 @@ struct Signals
    * have been included. Note that ConstraintMatrix::close() will be called
    * after this signal, so no need to call it inside.
    *
-   * The functions that can attach to this signal must take two
+   * By default, all the shared_ptrs points to the same
+   * ConstraintMatrix, which is constructed according to the boundary
+   * conditions set in the parameter file.
+   *
+   * The functions that can attach to this signal must take a
+   * std::vector<shared_ptr<ConstraintMatrix> > and a
    * ConstraintMatrix.
    */
-  boost::signals2::signal<void (ConstraintMatrix &constraints,
+  boost::signals2::signal<void (std::vector<shared_ptr<ConstraintMatrix> > &constraints,
                                 ConstraintMatrix &constraints_dot)> update_constraint_matrices;
 
   /**

source/pidomus.cc

@@ -58,6 +58,7 @@ piDoMUS<dim, spacedim, LAC>::piDoMUS (const std::string &name,
 
   pgr("Refinement"),
 
+  train_constraints(interface.n_matrices),
 
   current_time(std::numeric_limits<double>::quiet_NaN()),
   // IDA calls residual first with alpha = 0
@@ -107,6 +108,10 @@ piDoMUS<dim, spacedim, LAC>::piDoMUS (const std::string &name,
   interface.initialize_simulator (*this);
   lambdas.set_functions_to_default();
 
+  train_constraints[0] = SP(new ConstraintMatrix());
+
+  for (unsigned int i=1; i<n_matrices; ++i)
+    train_constraints[i] = train_constraints[0];
 
   for (unsigned int i=0; i<n_matrices; ++i)
     {
@@ -130,6 +135,7 @@ piDoMUS<dim, spacedim, LAC>::piDoMUS (const std::string &name,
 
   pgr("Refinement"),
 
+  train_constraints(interface.n_matrices),
 
   current_time(std::numeric_limits<double>::quiet_NaN()),
   // IDA calls residual first with alpha = 0
@@ -179,6 +185,7 @@ piDoMUS<dim, spacedim, LAC>::piDoMUS (const std::string &name,
   interface.initialize_simulator (*this);
   lambdas.set_functions_to_default();
 
+  train_constraints[0] = SP(new ConstraintMatrix);
 
   for (unsigned int i=0; i<n_matrices; ++i)
     {
@@ -207,7 +214,7 @@ void piDoMUS<dim, spacedim, LAC>::run ()
       else
         refine_mesh();
 
-      constraints.distribute(solution);
+      train_constraints[0]->distribute(solution);
       constraints_dot.distribute(solution_dot);
 
       if (time_stepper == "ida")
@@ -338,7 +345,7 @@ void piDoMUS<dim, spacedim, LAC>::setup_dofs (const bool &first_run)
       matrices[i]->clear();
       initializer(*matrix_sparsities[i],
                   *dof_handler,
-                  constraints,
+                  *train_constraints[i],
                   interface.get_matrix_coupling(i));
       matrices[i]->reinit(*matrix_sparsities[i]);
     }
@@ -353,8 +360,8 @@ void piDoMUS<dim, spacedim, LAC>::setup_dofs (const bool &first_run)
       else if (!we_are_parallel)
         {
           const QGauss<dim> quadrature_formula(fe->degree + 1);
-          VectorTools::project(interface.get_project_mapping(), *dof_handler, constraints, quadrature_formula, initial_solution, solution);
-          VectorTools::project(interface.get_project_mapping(), *dof_handler, constraints, quadrature_formula, initial_solution_dot, solution_dot);
+          VectorTools::project(interface.get_project_mapping(), *dof_handler, *train_constraints[0], quadrature_formula, initial_solution, solution);
+          VectorTools::project(interface.get_project_mapping(), *dof_handler, *train_constraints[0], quadrature_formula, initial_solution_dot, solution_dot);
         }
       else
         {

source/pidomus_assemble.cc

@@ -82,9 +82,9 @@ void piDoMUS<dim, spacedim, LAC>::assemble_matrices (const double t,
   {
 
     for (unsigned int i=0; i<n_matrices; ++i)
-      this->constraints.distribute_local_to_global (data.local_matrices[i],
-                                                    data.local_dof_indices,
-                                                    *(this->matrices[i]));
+      this->train_constraints[i]->distribute_local_to_global (data.local_matrices[i],
+                                                              data.local_dof_indices,
+                                                              *(this->matrices[i]));
   };
 
   auto local_assemble = [ this ]
@@ -217,9 +217,9 @@ piDoMUS<dim, spacedim, LAC>::residual(const double t,
   auto local_copy = [&dst, this] (const pidomus::CopyData & data)
   {
 
-    this->constraints.distribute_local_to_global (data.local_residual,
-                                                  data.local_dof_indices,
-                                                  dst);
+    this->train_constraints[0]->distribute_local_to_global (data.local_residual,
+                                                            data.local_dof_indices,
+                                                            dst);
   };
 
   auto local_assemble = [ this ]
@@ -257,7 +257,7 @@ piDoMUS<dim, spacedim, LAC>::residual(const double t,
   for (unsigned int i = 0; i < id.n_elements(); ++i)
     {
       auto j = id.nth_index_in_set(i);
-      if (constraints.is_constrained(j))
+      if (train_constraints[0]->is_constrained(j))
         dst[j] = solution(j) - locally_relevant_solution(j);
     }
 

source/pidomus_eigen.cc

@@ -78,7 +78,7 @@ void piDoMUS<dim, spacedim, LAC>::solve_eigenproblem()
       else
         refine_mesh();
 
-      constraints.distribute(solution);
+      train_constraints[0]->distribute(solution);
       constraints_dot.distribute(solution_dot);
 
       // initialize mass matrix
@@ -91,7 +91,7 @@ void piDoMUS<dim, spacedim, LAC>::solve_eigenproblem()
       mass_matrix.clear();
       initializer(mass_pattern,
                   *dof_handler,
-                  constraints,
+                  *train_constraints[0],
                   mass_coupling);
 
       mass_matrix.reinit(mass_pattern);
@@ -122,9 +122,9 @@ void piDoMUS<dim, spacedim, LAC>::solve_eigenproblem()
                         (const pidomus::CopyMass & data)
       {
 
-        this->constraints.distribute_local_to_global (data.local_matrix,
-                                                      data.local_dof_indices,
-                                                      mass_matrix);
+        this->train_constraints[0]->distribute_local_to_global (data.local_matrix,
+                                                                data.local_dof_indices,
+                                                                mass_matrix);
       };
 
       auto local_assemble = [ this ]

source/pidomus_helper_functions.cc

@@ -61,7 +61,7 @@ syncronize(const double &t,
       update_functions_and_constraints(t);
       typename LAC::VectorType tmp(solution);
       typename LAC::VectorType tmp_dot(solution_dot);
-      constraints.distribute(tmp);
+      train_constraints[0]->distribute(tmp);
       constraints_dot.distribute(tmp_dot);
 
       locally_relevant_solution = tmp;
@@ -84,7 +84,7 @@ syncronize(const double &t,
       update_functions_and_constraints(t);
       typename LAC::VectorType tmp(solution);
       typename LAC::VectorType tmp_dot(solution_dot);
-      constraints.distribute(tmp);
+      train_constraints[0]->distribute(tmp);
       constraints_dot.distribute(tmp_dot);
 
       locally_relevant_solution = tmp;
@@ -98,8 +98,9 @@ syncronize(const double &t,
 
       typename LAC::VectorType tmp(solution);
       typename LAC::VectorType tmp_dot(solution_dot);
-      constraints.distribute(tmp);
+      train_constraints[0]->distribute(tmp);
       constraints_dot.distribute(tmp_dot);
+
       locally_relevant_solution = tmp;
       locally_relevant_solution_dot = tmp_dot;
     }
@@ -113,7 +114,7 @@ syncronize(const double &t,
 
       typename LAC::VectorType tmp(solution);
       typename LAC::VectorType tmp_dot(solution_dot);
-      constraints.distribute(tmp);
+      train_constraints[0]->distribute(tmp);
       constraints_dot.distribute(tmp_dot);
 
       locally_relevant_solution = tmp;
@@ -134,31 +135,38 @@ void piDoMUS<dim, spacedim, LAC>::update_functions_and_constraints (const double
       forcing_terms.set_time(t);
       neumann_bcs.set_time(t);
     }
+
   // clear previously stored constraints
-  constraints.clear();
+  for (unsigned int i=0; i<n_matrices; ++i)
+    train_constraints[i]->clear();
   constraints_dot.clear();
 
   // compute hanging nodes
   DoFTools::make_hanging_node_constraints (*dof_handler,
-                                           constraints);
+                                           *train_constraints[0]);
 
   DoFTools::make_hanging_node_constraints (*dof_handler,
                                            constraints_dot);
 
-  // compute boundary values
-  apply_dirichlet_bcs(*dof_handler, dirichlet_bcs, constraints);
+  zero_average.apply_zero_average_constraints(*dof_handler, *train_constraints[0]);
+
+
+
+  // compute boundary values for the system matrix
+  apply_dirichlet_bcs(*dof_handler, dirichlet_bcs, *train_constraints[0]);
   apply_dirichlet_bcs(*dof_handler, dirichlet_bcs_dot, constraints_dot);
 
 
-  // apply zero average constraints
-  zero_average.apply_zero_average_constraints(*dof_handler, constraints);
+  // apply zero average constraints to the system matrix
+  zero_average.apply_zero_average_constraints(*dof_handler, *train_constraints[0]);
 
   // add user-supplied bcs
-  signals.update_constraint_matrices(constraints,constraints_dot);
+  signals.update_constraint_matrices(train_constraints, constraints_dot);
 
   // close the constraints
-  constraints.close ();
-  constraints_dot.close ();
+  for (unsigned int i=0; i<n_matrices; ++i)
+    train_constraints[i]->close();
+  constraints_dot.close();
 }
 
 template <int dim, int spacedim, typename LAC>
@@ -201,14 +209,14 @@ piDoMUS<dim, spacedim, LAC>::set_constrained_dofs_to_zero(typename LAC::VectorTy
   if (global_partitioning.n_elements() > 0)
     {
       auto k = global_partitioning.nth_index_in_set(0);
-      if (constraints.is_constrained(k))
+      if (train_constraints[0]->is_constrained(k))
         v[k] = 0;
       else
         v[k] = v[k];
       for (unsigned int i = 1; i < global_partitioning.n_elements(); ++i)
         {
           auto j = global_partitioning.nth_index_in_set(i);
-          if (constraints.is_constrained(j))
+          if (train_constraints[0]->is_constrained(j))
             v[j] = 0;
         }
       v.compress(VectorOperation::insert);

source/pidomus_mesh_refinement.cc

@@ -62,7 +62,7 @@ piDoMUS<dim, spacedim, LAC>::solver_should_restart(const double t,
       auto _timer = computing_timer.scoped_timer ("Compute error estimator");
       update_functions_and_constraints(t);
 
-      constraints.distribute(solution);
+      train_constraints[0]->distribute(solution);
       locally_relevant_solution = solution;
       constraints_dot.distribute(solution_dot);
       locally_relevant_solution_dot = solution_dot;
@@ -92,7 +92,7 @@ piDoMUS<dim, spacedim, LAC>::solver_should_restart(const double t,
                                         adaptive_refinement);
 
           update_functions_and_constraints(t);
-          constraints.distribute(solution);
+          train_constraints[0]->distribute(solution);
           constraints_dot.distribute(solution_dot);
 
           signals.fix_solutions_after_refinement(solution,solution_dot);
@@ -163,7 +163,7 @@ refine_and_transfer_solutions(LATrilinos::VectorType &y,
 
 
   update_functions_and_constraints(current_time);
-  constraints.distribute(y);
+  train_constraints[0]->distribute(y);
   constraints_dot.distribute(y_dot);
 
   locally_relevant_y = y;
@@ -213,10 +213,10 @@ refine_and_transfer_solutions(LADealII::VectorType &y,
   y_expl = new_sols[2];
 
   update_functions_and_constraints(previous_time);
-  constraints.distribute(y_expl);
+  train_constraints[0]->distribute(y_expl);
 
   update_functions_and_constraints(current_time);
-  constraints.distribute(y);
+  train_constraints[0]->distribute(y);
   constraints_dot.distribute(y_dot);
 
   locally_relevant_y = y;

tests/set_constrained_to_zero.cc

@@ -6,9 +6,9 @@ using namespace dealii;
 template<int fdim, int fspacedim, typename fn_LAC>
 void test(piDoMUS<fdim,fspacedim,fn_LAC>  &pi_foo)
 {
-  pi_foo.constraints.clear();
-  pi_foo.constraints.add_line(0);
-  pi_foo.constraints.close();
+  pi_foo.train_constraints[0]->clear();
+  pi_foo.train_constraints[0]->add_line(0);
+  pi_foo.train_constraints[0]->close();
 
   pi_foo.global_partitioning.clear();
   pi_foo.global_partitioning.set_size(2);