Dear Lucas,
thanks very much for replying.
I tried to solve the same problem without rotating the domain. I rotated
the problem.
Unfortunately the issue is still there.
I attach the new images.
I also attach the source code again because I realized that the previous
one had a misplaced bracket
and was giving an error at compile time. I attached the wrong version of
the file.
Thanks again.
Best,
Roberto
On Friday, April 13, 2018 at 7:54:54 AM UTC-4, Lucas Campos wrote:
>
> Dear Roberto,
>
> I think the problem is that you are using GridTools::transform (via
> GridTools::rotate), in a distributed triangulation. This should not be
> done, as noted here:
> https://www.dealii.org/8.5.0/doxygen/deal.II/namespaceGridTools.html#a212e99cf0d923cebfa04f1d23fa60b04
>
> If this is indeed the root of the issue, I see two options:
>
> 1. Rotate your problem, rather than your triangulation, 90 degrees.
> 2. You can create a temporary local normal triangulation, rotate it, and
> finally copy into the distributed triangulation. That is, something in the
> lines of
>
>
>>> Triangulation<2> tri_tmp;
>>
>> GridGenerator::hyper_cube_slit(tri_tmp, -1, 1, false);
>>
>> GridTools::rotate(-M_PI/2, tri_tmp);
>>
>>
>>> triangultion.copy_triangulation(tri_tmp);
>>
>> triangulation.refine_global (6);
>>
>> Bests,
> Lucas
> On Thursday, 12 April 2018 00:54:05 UTC+2, Roberto Porcù wrote:
>>
>> Dear all,
>>
>> in the reply to my first post I forgot the code.
>> I attach it to this reply.
>>
>> Best
>>
>> Roberto P.
>>
>>
>> On Tuesday, February 20, 2018 at 7:27:37 PM UTC-5, Roberto Porcù wrote:
>>>
>>> Dear all,
>>>
>>> I'm solving a linear elasticity problem on a cracked domain created by
>>> means
>>> of the GridGenerator::hyper_cube_slit function.
>>> When I solve it in serial everything works fine.
>>> When I solve it in parallel I get a problem on a node which is on one
>>> of the slit boundaries as it is possible to see from the attached
>>> pictures.
>>> I don't understand why that is happening because I wrote other parallel
>>> codes in deal.II
>>> and everything was always working fine. I've checked my code many times,
>>> also comparing It to the tutorials.
>>> I am supposing there is something related to the slit.
>>>
>>> Thanks in advance,
>>> Roberto
>>>
>>
--
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.
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/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.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/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.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_nothing.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/tria.h>
#include <fstream>
#include <iostream>
namespace CrackedDomain { namespace aux {
class AnalyticalSolution : public dealii::Function<2>
{
public:
AnalyticalSolution()
: dealii::Function<2>(1)
{}
virtual double
value(const dealii::Point<2> & p, const unsigned int component=0) const;
virtual dealii::Tensor<1,2>
gradient(const dealii::Point<2> & p, const unsigned int component=0) const;
virtual double
laplacian(const dealii::Point<2> & p, const unsigned int component=0) const;
virtual double
forcing(const dealii::Point<2> & p, const unsigned int component=0) const;
};
class LeftCrack : public dealii::Function<2>
{
public:
LeftCrack()
: dealii::Function<2>(1)
{}
virtual double
value(const dealii::Point<2> & p, const unsigned int component=0) const;
};
class RightCrack : public dealii::Function<2>
{
public:
RightCrack()
: dealii::Function<2>(1)
{}
virtual double
value(const dealii::Point<2> & p, const unsigned int component=0) const;
};
class LeftEdge : public dealii::Function<2>
{
public:
LeftEdge()
: dealii::Function<2>(1)
{}
virtual double
value(const dealii::Point<2> & p, const unsigned int component=0) const;
};
class RightEdge : public dealii::Function<2>
{
public:
RightEdge()
: dealii::Function<2>(1)
{}
virtual double
value(const dealii::Point<2> & p, const unsigned int component=0) const;
};
double
AnalyticalSolution::value(const dealii::Point<2> & p,
const unsigned int) const
{
const double x = p[0];
const double y = p[1];
if((y > 0) or (std::abs(y) < 1.e-14))
return 0;
if(x > 0)
return y*y*std::exp(-1*y);
if(x < 0)
return -1*y*y*std::exp(-1*y);
AssertThrow(false, dealii::ExcMessage("Error"));
return 0;
}
dealii::Tensor<1,2>
AnalyticalSolution::gradient(const dealii::Point<2> & p,
const unsigned int) const
{
dealii::Tensor<1,2> result;
result = 0;
const double x = p[0];
const double y = p[1];
if((y > 0) or (std::abs(y) < 1.e-14)) {
return result;
}
if(x > 0) {
result[1] = 2*y*std::exp(-1*y) - y*y*std::exp(-1*y);
return result;
}
if(x < 0) {
result[1] = -2*y*std::exp(-1*y) + y*y*std::exp(-1*y);
return result;
}
AssertThrow(false, dealii::ExcMessage("Error"));
return result;
}
double
AnalyticalSolution::laplacian(const dealii::Point<2> & p,
const unsigned int) const
{
const double x = p[0];
const double y = p[1];
if((y > 0) or (std::abs(y) < 1.e-14))
return 0;
if(x > 0)
return 2*std::exp(-1*y) - 4*y*std::exp(-1*y) + y*y*std::exp(-1*y);
if(x < 0)
return -2*std::exp(-1*y) + 4*y*std::exp(-1*y) - y*y*std::exp(-1*y);
AssertThrow(false, dealii::ExcMessage("Error"));
return 0;
}
double
AnalyticalSolution::forcing(const dealii::Point<2> & p,
const unsigned int component) const
{
return -1 * laplacian(p, component);
}
double
LeftCrack::value(const dealii::Point<2> & p,
const unsigned int) const
{
const double y = p[1];
return y*y*std::exp(-1*y);
}
double
RightCrack::value(const dealii::Point<2> & p,
const unsigned int) const
{
const double y = p[1];
return -1*y*y*std::exp(-1*y);
}
double
LeftEdge::value(const dealii::Point<2> &,
const unsigned int) const
{
return std::exp(1);
}
double
RightEdge::value(const dealii::Point<2> &,
const unsigned int) const
{
return -1*std::exp(1);
}
}}
namespace CrackedDomain
{
using namespace dealii;
class Problem
{
public:
Problem ();
~Problem ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void output_results () const;
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<2> triangulation;
DoFHandler<2> dof_handler;
FE_Q<2> fe;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
ConstraintMatrix constraints;
PETScWrappers::MPI::SparseMatrix system_matrix;
PETScWrappers::MPI::Vector locally_relevant_solution;
PETScWrappers::MPI::Vector system_rhs;
ConditionalOStream pcout;
};
Problem::Problem ()
: mpi_communicator (MPI_COMM_WORLD)
, triangulation (mpi_communicator)
, dof_handler (triangulation)
, fe (1)
, pcout (std::cout, Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
{}
Problem::~Problem ()
{
dof_handler.clear ();
}
void Problem::setup_system ()
{
dof_handler.distribute_dofs (fe);
locally_owned_dofs = dof_handler.locally_owned_dofs ();
DoFTools::extract_locally_relevant_dofs (dof_handler,
locally_relevant_dofs);
locally_relevant_solution.reinit (locally_owned_dofs,
locally_relevant_dofs, mpi_communicator);
system_rhs.reinit (locally_owned_dofs, mpi_communicator);
aux::AnalyticalSolution exact_sol;
aux::LeftCrack left_crack;
aux::RightCrack right_crack;
aux::LeftEdge left_edge;
aux::RightEdge right_edge;
constraints.clear ();
constraints.reinit (locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
for(unsigned int id(4); id < 7; ++id)
dealii::VectorTools::interpolate_boundary_values(
dof_handler, id, exact_sol, constraints);
dealii::VectorTools::interpolate_boundary_values(
dof_handler, 0, right_edge, constraints);
dealii::VectorTools::interpolate_boundary_values(
dof_handler, 1, left_crack, constraints);
dealii::VectorTools::interpolate_boundary_values(
dof_handler, 2, right_crack, constraints);
dealii::VectorTools::interpolate_boundary_values(
dof_handler, 3, left_edge, constraints);
constraints.close ();
DynamicSparsityPattern dsp (locally_relevant_dofs);
DoFTools::make_sparsity_pattern (dof_handler, dsp,
constraints, false);
SparsityTools::distribute_sparsity_pattern (dsp,
dof_handler.n_locally_owned_dofs_per_processor(),
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit (locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
void Problem::assemble_system ()
{
aux::AnalyticalSolution exact_sol;
const QGauss<2> quadrature_formula(2);
FEValues<2> 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> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
typename DoFHandler<2>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned()) {
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
const Point<2> & p = fe_values.quadrature_point(q_point);
const double rhs_value = exact_sol.forcing(p);
for (unsigned int i=0; i<dofs_per_cell; ++i) {
for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point) *
fe_values.JxW(q_point));
cell_rhs(i) += (rhs_value *
fe_values.shape_value(i,q_point) *
fe_values.JxW(q_point));
}
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
system_matrix.compress (VectorOperation::add);
system_rhs.compress (VectorOperation::add);
}
void Problem::solve ()
{
PETScWrappers::MPI::Vector
completely_distributed_solution (locally_owned_dofs, mpi_communicator);
SolverControl solver_control (dof_handler.n_dofs(), 1e-12);
PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
PETScWrappers::PreconditionBoomerAMG preconditioner;
PETScWrappers::PreconditionBoomerAMG::AdditionalData data;
data.symmetric_operator = true;
preconditioner.initialize(system_matrix, data);
solver.solve (system_matrix, completely_distributed_solution, system_rhs,
preconditioner);
pcout << " Solved in " << solver_control.last_step()
<< " iterations." << std::endl;
constraints.distribute (completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
}
void Problem::output_results () const
{
dealii::DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(locally_relevant_solution, "solution");
std::string filename;
const unsigned int n_processes =
dealii::Utilities::MPI::n_mpi_processes(mpi_communicator);
const unsigned int rank =
dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
if(n_processes > 1){
filename = "solution-rank"
+ dealii::Utilities::int_to_string(rank, 2)
+ ".vtu";
}
else{
filename = "solution.vtu";
}
data_out.build_patches();
std::ofstream output(filename);
data_out.write_vtu(output);
}
void Problem::run ()
{
GridGenerator::hyper_cube_slit(triangulation, -1, 1, true);
triangulation.refine_global (6);
// set boundary indicators
{
FE_Nothing<2> fe_nothing(2);
QMidpoint<1> quad;
FEFaceValues<2> fe_face_values(fe_nothing, quad, update_normal_vectors);
parallel::distributed::Triangulation<2>::active_cell_iterator
cell(triangulation.begin_active()),
endc(triangulation.end());
for(; cell != endc; ++cell) {
if(cell->is_locally_owned()) {
for(unsigned int f(0); f<GeometryInfo<2>::faces_per_cell; ++f) {
if(cell->face(f)->at_boundary()) {
fe_face_values.reinit(cell, f);
Point<2> face_center = cell->face(f)->center();
Tensor<1,2> n = fe_face_values.normal_vector(0);
const double linfty = std::max(std::abs(n[0]), std::abs(n[1]));
n /= linfty;
const double x = face_center[0];
const double y = face_center[1];
if((n[1] == -1) and (std::abs(y+1) < 1.e-13) and (x < 0))
cell->face(f)->set_all_boundary_ids(0);
if((n[0] == -1) and (std::abs(x) < 1.e-13))
cell->face(f)->set_all_boundary_ids(1);
if((n[0] == 1) and (std::abs(x) < 1.e-13))
cell->face(f)->set_all_boundary_ids(2);
if((n[1] == -1) and (std::abs(y+1) < 1.e-13) and (x > 0))
cell->face(f)->set_all_boundary_ids(3);
if((n[0] == 1) and (std::abs(x-1) < 1.e-13))
cell->face(f)->set_all_boundary_ids(4);
if((n[1] == 1) and (std::abs(y-1) < 1.e-13))
cell->face(f)->set_all_boundary_ids(5);
if((n[0] == -1) and (std::abs(x+1) < 1.e-13))
cell->face(f)->set_all_boundary_ids(6);
}
}
}
}
}
setup_system ();
pcout << " Number of active cells: "
<< triangulation.n_global_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
assemble_system ();
solve ();
if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32) {
output_results ();
}
pcout << std::endl;
}
}
int main(int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace CrackedDomain;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
Problem problem;
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;
}
