namespace internal
{
-#ifdef DEAL_II_WITH_P4EST
+#ifndef DEAL_II_WITH_P4EST
+ // Dummy implementation in case p4est is not available.
+ template <int dim,int spacedim>
+ class ExtrapolateImplementation
+ {
+ public:
+
+ ExtrapolateImplementation ()
+ {
+ Assert(false, ExcNotImplemented());
+ };
+
+ template <class InVector,class OutVector>
+ void extrapolate_parallel (const InVector &/*u2_relevant*/,
+ const DoFHandler<dim,spacedim> &/*dof2*/,
+ OutVector &/*u2*/)
+ {
+ Assert(false, ExcNotImplemented())
+ };
+ };
+#else
// Implementation of the @p extrapolate function
// on parallel distributed grids.
template <int dim,int spacedim>
namespace
{
+ template <class VectorType, class DH>
+ void reinit_distributed(const DH &/*dh*/,
+ VectorType &/*vector*/)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+#ifdef DEAL_II_WITH_PETSC
+ template <int dim, int spacedim>
+ void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+ PETScWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+ Assert (parallel_tria !=0, ExcNotImplemented());
+
+ const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+ vector.reinit(locally_owned_dofs, parallel_tria->get_communicator());
+ }
+#endif //DEAL_II_WITH_PETSC
+
+#ifdef DEAL_II_WITH_TRILINOS
+ template <int dim, int spacedim>
+ void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+ TrilinosWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+ Assert (parallel_tria !=0, ExcNotImplemented());
+
+ const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+ vector.reinit(locally_owned_dofs, parallel_tria->get_communicator());
+ }
+#endif //DEAL_II_WITH_TRILINOS
+
+ template <int dim, int spacedim, typename Number>
+ void reinit_distributed(const DoFHandler<dim, spacedim> &dh,
+ LinearAlgebra::distributed::Vector<Number> &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+ Assert (parallel_tria !=0, ExcNotImplemented());
+
+ const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+ vector.reinit(locally_owned_dofs, parallel_tria->get_communicator());
+ }
+
+
+
+ template <class VectorType, class DH>
+ void reinit_ghosted(const DH &/*dh*/,
+ VectorType &/*vector*/)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+#ifdef DEAL_II_WITH_PETSC
+ template <int dim, int spacedim>
+ void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
+ PETScWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+ Assert(parallel_tria !=0, ExcNotImplemented());
+ const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs);
+ vector.reinit (locally_owned_dofs, locally_relevant_dofs,
+ parallel_tria->get_communicator());
+ }
+#endif //DEAL_II_WITH_PETSC
+
+#ifdef DEAL_II_WITH_TRILINOS
+ template <int dim, int spacedim>
+ void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
+ TrilinosWrappers::MPI::Vector &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+ Assert (parallel_tria !=0, ExcNotImplemented());
+ const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs);
+ vector.reinit (locally_owned_dofs, locally_relevant_dofs,
+ parallel_tria->get_communicator());
+ }
+#endif //DEAL_II_WITH_TRILINOS
+
+ template <int dim, int spacedim, typename Number>
+ void reinit_ghosted(const DoFHandler<dim, spacedim> &dh,
+ LinearAlgebra::distributed::Vector<Number> &vector)
+ {
+ const parallel::distributed::Triangulation<dim,spacedim> *parallel_tria =
+ dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dh.get_tria());
+ Assert (parallel_tria !=0, ExcNotImplemented());
+ const IndexSet locally_owned_dofs = dh.locally_owned_dofs();
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs (dh, locally_relevant_dofs);
+ vector.reinit (locally_owned_dofs, locally_relevant_dofs,
+ parallel_tria->get_communicator());
+ }
+
+
+
template <int dim, class InVector, class OutVector, int spacedim>
void extrapolate_serial(const InVector &u3,
const DoFHandler<dim,spacedim> &dof2,
}
}
}
-
-
-
-#ifdef DEAL_II_WITH_P4EST
- template <int dim, class InVector, class OutVector, int spacedim>
- void extrapolate_parallel(const InVector &u2_relevant,
- const DoFHandler<dim,spacedim> &dof2,
- OutVector &u2)
- {
- if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dof2.get_tria()) != 0)
- {
- ExtrapolateImplementation<dim,spacedim> implementation;
- implementation.extrapolate_parallel (u2_relevant, dof2, u2);
- }
- else
- {
- extrapolate_serial (u2_relevant, dof2, u2);
- }
- }
-#endif //DEAL_II_WITH_P4EST
}
}
Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
- // make sure that each cell on the
- // coarsest level is at least once
- // refined, otherwise, these cells
- // can't be treated and would
- // generate a bogus result
+ // make sure that each cell on the coarsest level is at least once refined, otherwise, these cells
+ // can't be treated and would generate a bogus result
{
typename DoFHandler<dim,spacedim>::cell_iterator cell = dof2.begin(0),
endc = dof2.end(0);
for (; cell!=endc; ++cell)
- Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce());
+ Assert (cell->has_children() || cell->is_artificial(), ExcGridNotRefinedAtLeastOnce());
}
- OutVector u3;
- u3.reinit(u2);
- interpolate(dof1, u1, dof2, constraints, u3);
-
- internal::extrapolate_serial (u3, dof2, u2);
-
- constraints.distribute(u2);
- }
-
-
-#ifdef DEAL_II_WITH_P4EST
- template <int dim, class VectorType, int spacedim>
- void extrapolate_parallel(const DoFHandler<dim,spacedim> &dof1,
- const VectorType &u1,
- const DoFHandler<dim,spacedim> &dof2,
- const ConstraintMatrix &constraints2,
- VectorType &u2)
- {
- IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs();
- IndexSet dof2_locally_relevant_dofs;
- DoFTools::extract_locally_relevant_dofs (dof2,
- dof2_locally_relevant_dofs);
-
- VectorType u3 (dof2_locally_owned_dofs, u1.get_mpi_communicator ());
- interpolate (dof1, u1, dof2, constraints2, u3);
- VectorType u3_relevant (dof2_locally_owned_dofs,
- dof2_locally_relevant_dofs,
- u1.get_mpi_communicator ());
- u3_relevant = u3;
-
- internal::extrapolate_parallel (u3_relevant, dof2, u2);
-
- constraints2.distribute(u2);
- }
-
+ OutVector u3;
+ internal::reinit_distributed(dof2, u3);
+ if (dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim>*>(&dof2.get_tria()) != 0)
+ {
+ interpolate(dof1, u1, dof2, constraints, u3);
- template <int dim, int spacedim, typename Number>
- void extrapolate_parallel(const DoFHandler<dim,spacedim> &dof1,
- const LinearAlgebra::distributed::Vector<Number> &u1,
- const DoFHandler<dim,spacedim> &dof2,
- const ConstraintMatrix &constraints2,
- LinearAlgebra::distributed::Vector<Number> &u2)
- {
- IndexSet dof2_locally_owned_dofs = dof2.locally_owned_dofs();
- IndexSet dof2_locally_relevant_dofs;
- DoFTools::extract_locally_relevant_dofs (dof2,
- dof2_locally_relevant_dofs);
+ OutVector u3_relevant;
+ internal::reinit_ghosted(dof2, u3_relevant);
+ u3_relevant = u3;
- LinearAlgebra::distributed::Vector<Number> u3(dof2_locally_owned_dofs,
- dof2_locally_relevant_dofs,
- u1.get_mpi_communicator ());
- interpolate (dof1, u1, dof2, constraints2, u3);
- u3.update_ghost_values();
+ internal::ExtrapolateImplementation<dim,spacedim> implementation;
+ implementation.extrapolate_parallel (u3_relevant, dof2, u2);
+ }
+ else
+ {
+ interpolate(dof1, u1, dof2, constraints, u3);
- internal::extrapolate_parallel (u3, dof2, u2);
+ internal::extrapolate_serial (u3, dof2, u2);
+ }
- constraints2.distribute(u2);
+ constraints.distribute(u2);
}
-#endif //DEAL_II_WITH_P4EST
-
} // end of namespace FETools
#include <iomanip>
#include <string>
+//#define DEBUG_OUTPUT_VTK
template <int dim>
class TestFunction : public Function<dim>
const double global_error_before
= std::sqrt(Utilities::MPI::sum(std::pow(local_error_before,2.), MPI_COMM_WORLD));
- /*output_vector<dim, VectorType> (in_ghosted,
+#ifdef DEBUG_OUTPUT_VTK
+ output_vector<dim, VectorType> (in_ghosted,
Utilities::int_to_string(fe1.degree,1)+Utilities::int_to_string(dim,1)+std::string("in"),
- *dof1);*/
+ *dof1);
+#endif
{
const double l1_norm = in_distributed.l1_norm();
const double l2_norm = in_distributed.l2_norm();
<< linfty_norm << std::endl;
}
- FETools::extrapolate_parallel (*dof1, in_ghosted, *dof2, cm2, out_distributed);
+ FETools::extrapolate (*dof1, in_ghosted, *dof2, cm2, out_distributed);
out_distributed.compress(VectorOperation::insert);
out_ghosted = out_distributed;
-
- /*output_vector<dim, VectorType> (out_ghosted,
+#ifdef DEBUG_OUTPUT_VTK
+ output_vector<dim, VectorType> (out_ghosted,
Utilities::int_to_string(fe2.degree,1)+Utilities::int_to_string(dim,1)+std::string("out"),
- *dof2);*/
+ *dof2);
+#endif
{
const double l1_norm = out_distributed.l1_norm();
const double local_error_before = difference_before.l2_norm();
const double global_error_before
= std::sqrt(Utilities::MPI::sum(std::pow(local_error_before,2.), MPI_COMM_WORLD));
-
- /*output_vector<dim, VectorType> (in_ghosted,
+#ifdef DEBUG_OUTPUT_VTK
+ output_vector<dim, VectorType> (in_ghosted,
Utilities::int_to_string(fe1.degree,1)+Utilities::int_to_string(dim,1)+std::string("in"),
- *dof1);*/
+ *dof1);
+#endif
{
const double l1_norm = in_ghosted.l1_norm();
const double l2_norm = in_ghosted.l2_norm();
<< linfty_norm << std::endl;
}
- if (Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)==1)
- FETools::extrapolate (*dof1, in_ghosted, *dof2, cm2, out_ghosted);
- else
- FETools::extrapolate_parallel (*dof1, in_ghosted, *dof2, cm2, out_ghosted);
+ FETools::extrapolate (*dof1, in_ghosted, *dof2, cm2, out_ghosted);
out_ghosted.update_ghost_values();
-
- /*output_vector<dim, VectorType> (out_ghosted,
+#ifdef DEBUG_OUTPUT_VTK
+ output_vector<dim, VectorType> (out_ghosted,
Utilities::int_to_string(fe2.degree,1)+Utilities::int_to_string(dim,1)+std::string("out"),
- *dof2);*/
+ *dof2);
+#endif
{
const double l1_norm = out_ghosted.l1_norm();
const double l2_norm = out_ghosted.l2_norm();