I try to apply initial values to a vector defined as
*LinearAlgebraTrilinos::MPI::Vector
*using
VectorTools::project (dof_handler, hanging_node_constraints,
QGauss<dim>(fe.degree+1),
InitialValues<dim>(),
local_solution);
When initializing the variable fe (as FESystem<dim>) with one or two
components, it works fine. For more than two components I get the error
--------------------------------------------------------
An error occurred in line <1366> of file
<~/Downloads/dealii/include/deal.II/numerics/vector_tools.templates.h> in
function
void dealii::VectorTools::{anonymous}::project(const dealii::Mapping<dim>&,
const dealii::DoFHandler<dim>&, const dealii::ConstraintMatrix&, const
dealii::Quadrature<dim>&, const dealii::Function<dim, typename
VectorType::value_type>&, VectorType&, bool, const dealii
::Quadrature<(dim - 1)>&, bool) [with VectorType =
dealii::TrilinosWrappers::MPI::Vector; int dim = 2; typename
VectorType::value_type = double]
The violated condition was:
(dynamic_cast<const parallel::Triangulation<dim>* >
(&(dof.get_triangulation()))==nullptr)
Additional information:
You are trying to use functionality in deal.II that is currently not
implemented. In many cases, this indicates that there simply didn't appear much
of a need for it, or that the author of the original code did not have the time
to implement a particular case. If you
hit this exception, it is therefore worth the time to look into the code to
find out whether you may be able to implement the missing functionality. If you
do, please consider providing a patch to the deal.II development sources (see
the deal.II website on how to contri
bute).
Stacktrace:
-----------
#0 /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre:
#1 /opt/dealII/lib/libdeal_II.g.so.9.0.0-pre: void
dealii::VectorTools::project<2, dealii::TrilinosWrappers::MPI::Vector,
2>(dealii::DoFHandler<2, 2> const&, dealii::ConstraintMatrix const&,
dealii::Quadrature<2> const&, dealii::Function<2, dealii::TrilinosWrappers::MPI
::Vector::value_type> const&, dealii::TrilinosWrappers::MPI::Vector&, bool,
dealii::Quadrature<(2)-(1)> const&, bool)
#2 ./main: Step15::MinimalSurfaceProblem<2>::run()
#3 ./main: main
--------------------------------------------------------
[linux-lb8c:15830] *** Process received signal ***
[linux-lb8c:15830] Signal: Aborted (6)
[linux-lb8c:15830] Signal code: (-6)
[linux-lb8c:15830] [ 0] /lib64/libpthread.so.0(+0x12270)[0x7f294a477270]
[linux-lb8c:15830] [ 1] /lib64/libc.so.6(gsignal+0x110)[0x7f2946c1f0d0]
[linux-lb8c:15830] [ 2] /lib64/libc.so.6(abort+0x151)[0x7f2946c206b1]
[linux-lb8c:15830] [ 3]
/opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x6b9e5d1)[0x7f295b49e5d1]
[linux-lb8c:15830] [ 4]
/opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals5abortERKNS_13ExceptionBaseE+0x1a)[0x7f295b49edaf]
[linux-lb8c:15830] [ 5]
/opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii18deal_II_exceptions9internals11issue_errorINS_18StandardExceptions17ExcNotImplementedEEEvNS1_17ExceptionHandlingEPKciS7_S7_S7_T_+0x98)[0x7f2957373ea1]
[linux-lb8c:15830] [ 6]
/opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(+0x3f38e23)[0x7f2958838e23]
[linux-lb8c:15830] [ 7]
/opt/dealII/lib/libdeal_II.g.so.9.0.0-pre(_ZN6dealii11VectorTools7projectILi2ENS_16TrilinosWrappers3MPI6VectorELi2EEEvRKNS_10DoFHandlerIXT_EXT1_EEERKNS_16ConstraintMatrixERKNS_10QuadratureIXT_EEERKNS_8FunctionIXT1_ENT0_10value_typeEEERSH_bRKNSC_IX
miT_Li1EEEEb+0x2f)[0x7f295894906e]
[linux-lb8c:15830] [ 8]
./main(_ZN6Step1521MinimalSurfaceProblemILi2EE3runEv+0xc08)[0x420d08]
[linux-lb8c:15830] [ 9] ./main(main+0x3c)[0x414ad0]
[linux-lb8c:15830] [10]
/lib64/libc.so.6(__libc_start_main+0xea)[0x7f2946c09f4a]
[linux-lb8c:15830] [11] ./main(_start+0x2a)[0x41477a]
[linux-lb8c:15830] *** End of error message ***
Abgebrochen (Speicherabzug geschrieben)
when running in debug mode. It runs fine in release mode. Why does that happen
for more than two components, and how can I fix/circumvent that? Or did I
(again) forget something?
My minimal example is attached, the behaviour happens when setting
NUM_COMPONENTS via
#define NUM_COMPONENTS 100
to a value larger than 2.
Thank you!
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/data_out_base.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_sparsity_pattern.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <Sacado.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <experimental/type_traits>
#define NUM_COMPONENTS 100
template < typename T >
using toDim_t = decltype(std::declval<T>().norm());
template < typename T >
using has_toDim = std::experimental::is_detected< toDim_t, T >;
template <bool, typename T = void>
struct enable_if
{};
template <typename T>
struct enable_if<true, T> {
typedef T type;
};
namespace Step15
{
using namespace dealii;
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues() : Function<dim>(NUM_COMPONENTS * dim) {}
virtual double value(const Point<dim> &p, const unsigned int component) const;
virtual void vector_value(const Point<dim> &p, Vector<double> &values) const;
};
template <int dim>
double InitialValues<dim>::value(const Point<dim> &p, const unsigned int component) const
{
(void) p;
(void) component;
return 293;
}
template <int dim>
void InitialValues<dim>::vector_value(const Point<dim> &p, Vector<double> &values) const
{
for(size_t i = 0; i < values.size(); ++i)
values[i] = value(p, i);
}
template <int dim>
class MinimalSurfaceProblem
{
public:
MinimalSurfaceProblem ();
~MinimalSurfaceProblem ();
void run ();
private:
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
const size_t n_components;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
ConstraintMatrix hanging_node_constraints;
LinearAlgebraTrilinos::MPI::Vector present_solution;
LinearAlgebraTrilinos::MPI::Vector old_solution;
};
template <int dim>
MinimalSurfaceProblem<dim>::MinimalSurfaceProblem ()
:
mpi_communicator(MPI_COMM_WORLD),
triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(Triangulation<dim>::smoothing_on_refinement | Triangulation<dim>::smoothing_on_coarsening)),
n_components(NUM_COMPONENTS),
dof_handler (triangulation),
fe (FE_Q<dim>(2), n_components * dim)
{
}
template <int dim>
MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem ()
{
dof_handler.clear ();
}
template <int dim>
void MinimalSurfaceProblem<dim>::run ()
{
GridGenerator::hyper_cube(triangulation, -1., 1.);
for(auto cell = triangulation.begin(); cell != triangulation.end(); ++cell)
for(size_t face_number = 0; face_number < GeometryInfo<dim>::faces_per_cell; ++face_number)
cell->face(face_number)->set_boundary_id(0);
triangulation.refine_global(3);
dof_handler.distribute_dofs (fe);
const size_t dof_numbers = dof_handler.n_dofs();
IndexSet solution_partitioning(dof_numbers), solution_relevant_partitioning(dof_numbers);
solution_partitioning = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
hanging_node_constraints.clear();
hanging_node_constraints.reinit(solution_relevant_partitioning);
DoFTools::make_hanging_node_constraints(dof_handler, hanging_node_constraints);
hanging_node_constraints.close();
LinearAlgebraTrilinos::MPI::Vector local_solution;
local_solution.reinit(dof_handler.locally_owned_dofs(), mpi_communicator);
DoFTools::extract_locally_relevant_dofs(dof_handler, solution_relevant_partitioning);
VectorTools::project (dof_handler,
hanging_node_constraints,
QGauss<dim>(fe.degree+1),
InitialValues<dim>(),
local_solution);
//boundary_constraints.distribute(local_solution);
hanging_node_constraints.distribute(local_solution);
old_solution = local_solution;
present_solution = local_solution;
}
}
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step15;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
MinimalSurfaceProblem<2> laplace_problem_2d;
laplace_problem_2d.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
