Hi all,
I have a problem when using GridGenerator::cylinder_shell. I have a mesh
and when I refine it I want to have a cylinder but when the code runs I
have a seg fault with GridGenerator::cylinder_shell :(
Can you tell me what is wrong in my code please ?
Here is the function mesh_refinement and I join files Mesh.geo and
laplacian.cc
void LaplaceProblem::mesh_refinement(){
int Rint=241, Rext=240 ; //boundaries indicators
double inner_radius_value=0.2;
double outer_radius_value=0.3;
double length = 0.04;
GridGenerator::cylinder_shell
(triangulation,length,inner_radius_value,outer_radius_value);
int x=0,y=1,z=2;
static const CylinderBoundary<3> inner_cylinder (inner_radius_value,
z);
static const CylinderBoundary<3> outer_cylinder (outer_radius_value,
z);
triangulation.set_boundary (Rint, inner_cylinder);
triangulation.set_boundary (Rext, outer_cylinder);
triangulation.refine_global(1);
}
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/grid_in.h>
#include <grid/grid_reordering.h>
#include <grid/grid_tools.h>
#include <grid/grid_refinement.h>
#include <grid/tria_boundary_lib.h>
#include <fe/fe_q.h>
#include <fe/fe_tools.h>
#include <fe/fe_values.h>
#include <dofs/dof_tools.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <base/quadrature_lib.h>
#include <base/function.h>
#include <base/parameter_handler.h>
#include <base/path_search.h>
#include <base/parsed_function.h>
#include <base/utilities.h>
#include <lac/vector.h>
#include <lac/full_matrix.h>
#include <lac/sparse_matrix.h>
#include <lac/solver_cg.h>
#include <lac/precondition.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_vector.h>
#include <lac/trilinos_precondition.h>
#include <lac/trilinos_solver.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_block_vector.h>
#include <lac/trilinos_precondition.h>
#include <lac/full_matrix.h>
#include <lac/solver_gmres.h>
#include <lac/solver_cg.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/constraint_matrix.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <numeric>
#include <fstream>
#include <iostream>
#include <fstream>
#include <string>
#include <ctype.h>
#include <stdlib.h>
#include <map>
#include <memory>
using namespace dealii;
class LaplaceProblem {
public:
LaplaceProblem ();
void run ();
private:
void initialization ();
void assemble_system ();
void solve ();
void output_results () const;
void Kelly_estimate_error(Vector<float> &error_indicators);
void mesh_refinement();
void read_mesh();
Triangulation<3> triangulation;
FE_Q<3> fe;
DoFHandler<3> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> Matrix;
Vector<double> solution;
Vector<double> Rhs;
};
LaplaceProblem::LaplaceProblem ()
:
fe (1), // Q1
dof_handler (triangulation)
{}
void LaplaceProblem::read_mesh ()
{
PathSearch search("./");
std::string file;
file = search.find("Mesh.msh");
std::ifstream fichier(file.c_str());
GridIn<3> grid;
grid.attach_triangulation(triangulation);
grid.read_msh(fichier);
//triangulation.refine_global (7);
}
void LaplaceProblem::initialization ()
{
dof_handler.distribute_dofs (fe);
sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
sparsity_pattern.compress ();
Matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
Rhs.reinit (dof_handler.n_dofs());
}
//solve a simple problem :
// -Delta u = 0
// u = 0.25 on top
// u = 0 on bottom
// \nabla u \cdot n = 0 on other boundaries
// solution is given by u= \tehta / (2 pi )
void LaplaceProblem::assemble_system ()
{
int nb_quadrature_points = 2;
QGauss<3> quadrature_formula(nb_quadrature_points);
QGauss<2> face_quadrature_formula(nb_quadrature_points);
FEValues<3> fe_values (fe, quadrature_formula,
update_values |
update_gradients |
update_hessians |
update_JxW_values|
update_q_points);
FEFaceValues<3> ffv (fe, face_quadrature_formula,
update_values |
update_gradients |
update_JxW_values |
update_normal_vectors |
update_q_points);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int Nq = quadrature_formula.size();
const unsigned int Nqf = face_quadrature_formula.size();
int bottom=201, top=101, Rint=241, Rext=240 , channel=52; //boundaries indicators
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
DoFHandler<3>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
cell->get_dof_indices (local_dof_indices);
for(unsigned int q=0; q<Nq; q++){
for (unsigned int i=0; i<dofs_per_cell; ++i){
for (unsigned int j=0; j<dofs_per_cell; ++j){
Matrix.add (local_dof_indices[i],
local_dof_indices[j],
fe_values.shape_grad(i,q) *
fe_values.shape_grad(j,q) *
fe_values.JxW (q));
}//end j
}//end i
}//end q
double gamma=100; // penaldir
Point<3> N;
// Dirichlet in a weak form
for (unsigned int face_no=0; face_no<GeometryInfo<3>::faces_per_cell; ++face_no)
{
int face_indicator = cell->face(face_no)->boundary_indicator();
if(face_indicator==top || face_indicator==bottom )
{
ffv.reinit (cell, face_no);
const std::vector<Point<3> > &q_points = ffv.get_quadrature_points();
for (unsigned int q=0; q<Nqf; ++q)
{
double v = 0.25*(face_indicator==top);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
Matrix.add (local_dof_indices[i],
local_dof_indices[j],
-ffv.shape_grad(i,q)*N* ffv.shape_value(j,q) * ffv.JxW(q)
-ffv.shape_grad(j,q)*N* ffv.shape_value(i,q) * ffv.JxW(q)
+ gamma/cell->face(face_no)->diameter() *ffv.shape_value(i,q)*ffv.shape_value(j,q) *ffv.JxW(q) );
Rhs(local_dof_indices[i]) += (v*gamma/cell->face(face_no)->diameter())*ffv.shape_value(i,q)*ffv.JxW(q)- ffv.shape_grad(i,q)*N*v *ffv.JxW(q);
} // end of loop over i
} //end of loop over quadrature points
} //end of the test
} //end of loop over faces
} //end of loop over cells
}
void LaplaceProblem::mesh_refinement(){
int Rint=241, Rext=240 ; //boundaries indicators
double inner_radius_value=0.2;
double outer_radius_value=0.3;
double length = 0.04;
GridGenerator::cylinder_shell (triangulation,length,inner_radius_value,outer_radius_value);
int x=0,y=1,z=2;
static const CylinderBoundary<3> inner_cylinder (inner_radius_value, z);
static const CylinderBoundary<3> outer_cylinder (outer_radius_value, z);
triangulation.set_boundary (Rint, inner_cylinder);
triangulation.set_boundary (Rext, outer_cylinder);
triangulation.refine_global(1);
}
void LaplaceProblem::solve ()
{
SolverControl solver_control (10000, 1e-8);
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize (Matrix,1.2);
cg.solve (Matrix, solution, Rhs, preconditioner);
}
void LaplaceProblem::output_results(void) const
{
DataOut<3> data_out_solution;
data_out_solution.attach_dof_handler (dof_handler);
data_out_solution.add_data_vector (solution, "solution");
data_out_solution.build_patches ();
std::ostringstream filename;
filename << "solution.vtk";
std::ofstream output (filename.str().c_str());
data_out_solution.write_vtk (output);
}
void LaplaceProblem::run ()
{
read_mesh();
initialization();
assemble_system ();
solve ();
output_results() ;
mesh_refinement();
initialization ();
assemble_system ();
solve ();
}
int main ()
{
LaplaceProblem laplace_problem;
laplace_problem.run ();
return 0;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
