Dear all,
I am trying to solve a poroelastic problem by using a FE_System consisting
of FE_Q for displacements, FE_RaviartThomas for velocities and FE_DGQ for
pressure and I want to use an academic example to check that the code is
correct. This example needs to implement a non-null right-hand-side vector.
And here it is my problem, when compiling my code, I obtain an error when
defining the local right-hand-side vector.
Here you can see the make error I have obtained:
/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc: In instantiation of
‘*void Poroelastic::PoroelasticProblem<dim>::assemble_system() [with int
dim = 2]*’:
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:593:21:* required
from ‘*void Poroelastic::PoroelasticProblem<dim>::run() [with int dim = 2]*’
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:612:31:* required
from here
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:545:37:* *error: *no
match for ‘*operator**’ (operand types are
‘*__gnu_cxx::__alloc_traits<std::allocator<dealii::Vector<double>
> >::value_type {aka dealii::Vector<double>}*’ and ‘*const
dealii::Tensor<1, 2>*’)
(disp_rhs_values[q] * phi_i_u - pres_rhs_values[q] *
phi_i_p) * fe_values.JxW(q);
* ^*
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:545:37:* *note:
*candidates
are:
In file included from
*/home/teresa/.local/bin/deal.II/include/deal.II/lac/block_vector_base.h:29:0*
,
from
*/home/teresa/.local/bin/deal.II/include/deal.II/lac/block_vector.h:25*,
from
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:30*:
*/home/teresa/.local/bin/deal.II/include/deal.II/lac/vector.h:478:10:*
*note: *template<class Number2> Number
dealii::Vector<Number>::operator*(const dealii::Vector<OtherNumber>&) const
[with Number2 = Number2; Number = double]
Number operator*(const Vector<Number2> &V) const;
* ^*
*/home/teresa/.local/bin/deal.II/include/deal.II/lac/vector.h:478:10:*
*note: * template argument deduction/substitution failed:
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:545:37:* *note: *
‘*const
dealii::Tensor<1, 2>*’ is not derived from ‘*const
dealii::Vector<OtherNumber>*’
(disp_rhs_values[q] * phi_i_u - pres_rhs_values[q] *
phi_i_p) * fe_values.JxW(q);
* ^*
In file included from
*/home/teresa/.local/bin/deal.II/include/deal.II/base/template_constraints.h:22:0*
,
from
*/home/teresa/.local/bin/deal.II/include/deal.II/base/tensor.h:24*,
from
*/home/teresa/.local/bin/deal.II/include/deal.II/base/point.h:23*,
from
*/home/teresa/.local/bin/deal.II/include/deal.II/base/quadrature.h:22*,
from
*/home/teresa/.local/bin/deal.II/include/deal.II/base/quadrature_lib.h:22*,
from
*/home/teresa/dealii/Poroelast_simple/Poroelasticity.cc:25*:
*/home/teresa/.local/bin/deal.II/include/deal.II/base/complex_overloads.h:41:1:*
*note: *template<class T, class U> typename
dealii::ProductType<std::complex<_Tp>, std::complex<_Up> >::type
dealii::operator*(const std::complex<_Tp>&, const std::complex<_Up>&)
operator*(const std::complex<T> &left, const std::complex<U> &right)
* ^*
Maybe the make error is a key to find the error, but I am unable to
understand the problem.
I send attached a simplified version of the code that I am using.
Does anyone know what is the problem?
Thank you in advance.
Best regards,
Teresa.
--
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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4f229649-7775-4bff-ab96-b7554dcbdbd7%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2006 - 2019 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
*/
// This program is an adaptation of step-21 to solve a poroelastic problem
// @sect3{Include files}
// All of these include files have been used before:
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparse_direct.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/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <iostream>
#include <fstream>
#include <deal.II/base/tensor_function.h>
namespace Poroelastic
{
using namespace dealii;
template <int dim>
class PoroelasticProblem
{
public:
PoroelasticProblem(const unsigned int degree);
void run();
private:
void make_grid_and_dofs();
void assemble_system();
void solve();
void output_results() const;
const unsigned int degree;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
ConstraintMatrix constraints;
const unsigned int n_refinement_steps;
double time_step;
unsigned int timestep_number;
double end_time;
double lambda;
double mu;
double present_time;
double viscosity;
BlockVector<double> solution;
BlockVector<double> system_rhs;
};
// @sect3{Equation data}
// @sect4{Pressure right hand side}
template <int dim>
class PressureRightHandSide : public Function<dim>
{
public:
PressureRightHandSide()
: Function<dim>(1)
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
};
template <int dim>
double PressureRightHandSide<dim>::value(const Point<dim> &p,
const unsigned int) const
{
double present_time=1.0;
return 2.*(4.*numbers::PI*std::sin(2*numbers::PI*present_time)+
std::cos(2*numbers::PI*present_time))*numbers::PI*std::sin(2*numbers::PI*p[0])*std::sin(2*numbers::PI*p[1]);
}
// @sect4{Exact Solution}
template <int dim>
class ExactSolution : public Function <dim>
{
public:
ExactSolution(const double present_time);
virtual void vector_value(const Point<dim> &p,
Vector<double> & values) const override;
virtual void vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const override;
private:
const double present_time;
};
template <int dim>
ExactSolution<dim>::ExactSolution(const double present_time)
: Function<dim>(dim + dim +1)
, present_time(present_time)
{}
template <int dim>
void ExactSolution<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
Assert(values.size() == dim + 1,
ExcDimensionMismatch(values.size(), dim + 1));
values(0) = -std::cos(2*numbers::PI*p[0])*std::sin(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)/(4*numbers::PI);
values(1) = -std::sin(2*numbers::PI*p[0])*std::cos(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)/(4*numbers::PI);
values(2) = -2.0*std::cos(2*numbers::PI*p[0])*std::sin(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)*numbers::PI;
values(3) = -2.0*std::sin(2*numbers::PI*p[0])*std::cos(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)*numbers::PI;
values(4) = -2.0*std::sin(2*numbers::PI*p[0])*std::cos(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)*numbers::PI;
}
template <int dim>
void ExactSolution<dim>::vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const
{
const unsigned int n_points = points.size();
Assert(value_list.size() == n_points,ExcDimensionMismatch(value_list.size(), n_points));
for (unsigned int p = 0; p < n_points; ++p)
ExactSolution<dim>::vector_value(points[p], value_list[p]);
}
// @sect4{Displacements right hand side}
template <int dim>
class DispRightHandSide: public Function <dim>
{
public:
//DispRightHandSide(const double present_time);
DispRightHandSide();
virtual void vector_value(const Point<dim> &p,
Vector<double> & values) const override;
virtual void vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const override;
//private:
// const double present_time;
};
template <int dim>
//DispRightHandSide<dim>::DispRightHandSide(const double present_time)
DispRightHandSide<dim>::DispRightHandSide()
: Function<dim>(dim)
//, present_time(present_time)
{}
template <int dim>
void DispRightHandSide<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
Assert(values.size() == dim + 1,
ExcDimensionMismatch(values.size(), dim + 1));
double present_time = 1.0;
values = 0;
values(0) = -4.*numbers::PI*std::sin(2*numbers::PI*present_time)*std::cos(2*numbers::PI*p[0])*std::sin(2*numbers::PI*p[1]);
values(1) = -4.*numbers::PI*std::sin(2*numbers::PI*present_time)*std::sin(2*numbers::PI*p[0])*std::cos(2*numbers::PI*p[1]);
}
template <int dim>
void DispRightHandSide<dim>::vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const
{
const unsigned int n_points = points.size();
Assert(value_list.size() == n_points,ExcDimensionMismatch(value_list.size(), n_points));
for (unsigned int p = 0; p < n_points; ++p)
DispRightHandSide<dim>::vector_value(points[p], value_list[p]);
}
// @sect4{Displacement Boundary Values}
//
template <int dim>
class DisplacementBoundaryValues : public Function<dim>
{
public:
DisplacementBoundaryValues(const double present_time);
virtual void vector_value(const Point<dim> &p,
Vector<double> & values) const override;
virtual void vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const override;
private:
const double present_time;
};
template <int dim>
DisplacementBoundaryValues<dim>::DisplacementBoundaryValues(const double present_time)
: Function<dim>(dim + dim +1)
, present_time(present_time)
{}
template <int dim>
void DisplacementBoundaryValues<dim>::vector_value(const Point<dim> & p,
Vector<double> &values) const
{
Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
values = 0;
values(0) = -std::cos(2*numbers::PI*p[0])*std::sin(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)/(4*numbers::PI);
values(1) = -std::sin(2*numbers::PI*p[0])*std::cos(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)/(4*numbers::PI);
}
template <int dim>
void DisplacementBoundaryValues<dim>::vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const
{
const unsigned int n_points = points.size();
Assert(value_list.size() == n_points,ExcDimensionMismatch(value_list.size(), n_points));
for (unsigned int p = 0; p < n_points; ++p)
DisplacementBoundaryValues<dim>::vector_value(points[p], value_list[p]);
}
// @sect4{Velocity Boundary Values}
//
template <int dim>
class VelocityBoundaryValues : public Function<dim>
{
public:
VelocityBoundaryValues(const double present_time);
virtual void vector_value(const Point<dim> &p,
Vector<double> & values) const override;
virtual void vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const override;
private:
const double present_time;
};
template <int dim>
VelocityBoundaryValues<dim>::VelocityBoundaryValues(const double present_time)
: Function<dim>(dim + dim +1)
, present_time(present_time)
{}
template <int dim>
void VelocityBoundaryValues<dim>::vector_value(const Point<dim> & p,
Vector<double> &values) const
{
Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
values = 0;
values(2) = -2.0*std::cos(2*numbers::PI*p[0])*std::sin(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)*numbers::PI;
values(3) = -2.0*std::sin(2*numbers::PI*p[0])*std::cos(2*numbers::PI*p[1])*std::sin(2*numbers::PI*present_time)*numbers::PI;
}
template <int dim>
void VelocityBoundaryValues<dim>::vector_value_list(const std::vector<Point<dim>> &points,
std::vector<Vector<double>> & value_list) const
{
const unsigned int n_points = points.size();
Assert(value_list.size() == n_points,ExcDimensionMismatch(value_list.size(), n_points));
for (unsigned int p = 0; p < n_points; ++p)
VelocityBoundaryValues<dim>::vector_value(points[p], value_list[p]);
}
// @sect4{Initial data}
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues()
: Function<dim>(dim + dim + 1)
{}
virtual double value(const Point<dim> & p,
const unsigned int component = 0) const override
{
return Functions::ZeroFunction<dim>(dim + dim + 1).value(p, component);
}
virtual void vector_value(const Point<dim> &p,
Vector<double> & values) const override
{
Functions::ZeroFunction<dim>(dim + dim + 1).vector_value(p, values);
}
};
// @sect3{The inverse permeability tensor}
namespace IdentityPermeability
{
template <int dim>
class KInverse : public TensorFunction<2, dim>
{
public:
KInverse()
: TensorFunction<2, dim>()
{}
virtual void
value_list(const std::vector<Point<dim>> &points,
std::vector<Tensor<2, dim>> & values) const override
{
Assert(points.size() == values.size(),
ExcDimensionMismatch(points.size(), values.size()));
for (unsigned int p = 0; p < points.size(); ++p)
{
values[p].clear();
const double permeability = 1.0;
for (unsigned int d = 0; d < dim; ++d)
values[p][d][d] = 1./permeability;
}
}
};
}
// @sect3{<code>PoroelasticProblem</code> class implementation}
// @sect4{PoroelasticProblem::PoroelasticProblem}
template <int dim>
PoroelasticProblem<dim>::PoroelasticProblem(const unsigned int degree)
: degree(degree)
, fe(FE_Q<dim>(degree+1),dim,
FE_RaviartThomas<dim>(degree),
1,
FE_DGQ<dim>(degree),
1)
, dof_handler(triangulation)
//, n_refinement_steps(5)
, n_refinement_steps(1)
, time_step(1)
, timestep_number(1)
, end_time(1.0)
, lambda(1)
, mu(1)
, present_time(0.0)
, viscosity(0.2)
{}
// @sect4{PoroelasticProblem::make_grid_and_dofs}
template <int dim>
void PoroelasticProblem<dim>::make_grid_and_dofs()
{
GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(n_refinement_steps);
dof_handler.distribute_dofs(fe);
DoFRenumbering::component_wise(dof_handler);
std::vector<types::global_dof_index> dofs_per_component(dim + dim + 1);
DoFTools::count_dofs_per_component(dof_handler, dofs_per_component);
const unsigned int n_u = dofs_per_component[0]+dofs_per_component[1],
n_v = dofs_per_component[dim],
n_p = dofs_per_component[dim + dim];
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< "Number of active vertices: " << triangulation.n_used_vertices()
<< std::endl
<< "Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_u << '+' << n_v << '+' << n_p << ')' << std::endl
<< std::endl;
const unsigned int n_couplings = dof_handler.max_couplings_between_dofs();
sparsity_pattern.reinit(3, 3);
sparsity_pattern.block(0, 0).reinit(n_u, n_u, n_couplings);
sparsity_pattern.block(1, 0).reinit(n_v, n_u, n_couplings);
sparsity_pattern.block(2, 0).reinit(n_p, n_u, n_couplings);
sparsity_pattern.block(0, 1).reinit(n_u, n_v, n_couplings);
sparsity_pattern.block(1, 1).reinit(n_v, n_v, n_couplings);
sparsity_pattern.block(2, 1).reinit(n_p, n_v, n_couplings);
sparsity_pattern.block(0, 2).reinit(n_u, n_p, n_couplings);
sparsity_pattern.block(1, 2).reinit(n_v, n_p, n_couplings);
sparsity_pattern.block(2, 2).reinit(n_p, n_p, n_couplings);
sparsity_pattern.collect_sizes();
DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
sparsity_pattern.compress();
system_matrix.reinit(sparsity_pattern);
solution.reinit(3);
solution.block(0).reinit(n_u);
solution.block(1).reinit(n_v);
solution.block(2).reinit(n_p);
solution.collect_sizes();
system_rhs.reinit(3);
system_rhs.block(0).reinit(n_u);
system_rhs.block(1).reinit(n_v);
system_rhs.block(2).reinit(n_p);
system_rhs.collect_sizes();
}
// @sect4{PoroelasticProblem::assemble_system}
template <int dim>
void PoroelasticProblem<dim>::assemble_system()
{
system_matrix = 0;
system_rhs = 0;
QGauss<dim> quadrature_formula(degree + 2);
QGauss<dim - 1> face_quadrature_formula(degree + 2);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const DispRightHandSide<dim> disp_right_hand_side;
const PressureRightHandSide<dim> pres_right_hand_side;
const IdentityPermeability::KInverse<dim> k_inverse;
std::vector<Vector<double>> disp_rhs_values(n_q_points,Vector<double>(dim));
std::vector<double> pres_rhs_values(n_q_points);
std::vector<Tensor<2, dim>> k_inverse_values(n_q_points);
const FEValuesExtractors::Vector displacements(0);
const FEValuesExtractors::Vector velocities(dim);
const FEValuesExtractors::Scalar pressure(2*dim);
ComponentMask disp_mask = fe.component_mask(displacements);
ComponentMask vel_mask = fe.component_mask(velocities);
std::vector<Tensor<2, dim>> elasticity_grad_phi(dofs_per_cell);
std::vector<double> elasticity_div_phi(dofs_per_cell);
std::vector<Tensor<1, dim>> elasticity_phi(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
local_matrix = 0;
local_rhs = 0;
disp_right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
disp_rhs_values);
pres_right_hand_side.value_list(fe_values.get_quadrature_points(),
pres_rhs_values);
k_inverse.value_list(fe_values.get_quadrature_points(),
k_inverse_values);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
elasticity_grad_phi[k] =
fe_values[displacements].gradient(k, q);
elasticity_div_phi[k] =
fe_values[displacements].divergence(k, q);
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const Tensor<1, dim> phi_i_u = fe_values[displacements].value(i, q);
const Tensor<1, dim> phi_i_v = fe_values[velocities].value(i, q);
const double div_phi_i_u = fe_values[displacements].divergence(i, q);
const double div_phi_i_v = fe_values[velocities].divergence(i, q);
const double phi_i_p = fe_values[pressure].value(i, q);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const Tensor<1, dim> phi_j_v = fe_values[velocities].value(j, q);
const double div_phi_j_u = fe_values[displacements].divergence(j, q);
const double div_phi_j_v = fe_values[velocities].divergence(j, q);
const double phi_j_p = fe_values[pressure].value(j, q);
local_matrix(i, j) +=
(lambda * elasticity_div_phi[i] * elasticity_div_phi[j] +
mu * scalar_product(elasticity_grad_phi[i],
elasticity_grad_phi[j]) +
mu * scalar_product(elasticity_grad_phi[i],
transpose(elasticity_grad_phi[j]))
-div_phi_i_u * phi_j_p +
phi_i_v * k_inverse_values[q] * phi_j_v-
div_phi_i_v * phi_j_p
-div_phi_j_u * phi_i_p - div_phi_j_v * phi_i_p
)*fe_values.JxW(q);
}
local_rhs(i) +=
(disp_rhs_values[q] * phi_i_u - pres_rhs_values[q] * phi_i_p) * fe_values.JxW(q);
}
}
// The final step in the loop over all cells is to transfer local
// contributions into the global matrix and right hand side vector:
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
local_matrix(i, j));
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
}
std::map<types::global_dof_index, double> disp_boundary_values;
std::map<types::global_dof_index, double> vel_boundary_values;
//constraints.clear();
VectorTools::interpolate_boundary_values(dof_handler,
0,
DisplacementBoundaryValues<dim>(present_time),
disp_boundary_values,disp_mask);
VectorTools::interpolate_boundary_values(dof_handler,
0,
VelocityBoundaryValues<dim>(present_time),
vel_boundary_values,vel_mask);
//constraints.close();
MatrixTools::apply_boundary_values(disp_boundary_values,
system_matrix,
solution,
system_rhs);
MatrixTools::apply_boundary_values(vel_boundary_values,
system_matrix,
solution,
system_rhs);
}
// @sect4{PoroelasticProblem::solve}
// @sect4{TwoPhaseFlowProblem::run}
template <int dim>
void PoroelasticProblem<dim>::run()
{
make_grid_and_dofs();
assemble_system();
}
} // namespace Poroelastic
// @sect3{The <code>main</code> function}
// That's it. In the main function, we pass the degree of the finite element
// space to the constructor of the PoroelasticProblem object. Here, we use
// zero-th degree elements, i.e. $RT_0\times DQ_0 \times DQ_0$. The rest is as
// in all the other programs.
int main()
{
try
{
using namespace dealii;
using namespace Poroelastic;
PoroelasticProblem<2> poroelastic_problem(0);
poroelastic_problem.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;
}
