Sorry, I have sent the previous mail by mistake... Now I join files
correctly (laplacian.cc is a short 3D problem)
Thank you in advance
Bests,
Stéphane Veys
#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;
}
/******************************************************************/
/* Geometrie pour 2D_Modelization_longitudinally_cooled_helices */
/******************************************************************/
/* nombre de points suivant r et theta */
nr = 5;
ntheta = 20;
nz = 2;
//nr = 50;
//ntheta = 50;
/* rayons externes et internes */
r1=0.2;
r2=0.3;
Lc=1e-2;
Point(1001) = {0, 0 , 0.0 , Lc};
Point(1002) = { r1, 0.0 , 0.0 , Lc};
Point(1003) = { r2, 0.0 , 0.0 , Lc};
Point(1004) = { 0, r2 , 0.0 , Lc};
Point(1005) = { 0, r1 , 0.0 , Lc};
Circle(1007) = {1003,1001,1004};
Circle(1008) = {1005,1001,1002};
Line(1009) = {1004,1005};
Line(1010) = {1002,1003};
Line Loop(1011) = {1010,1007,1009,1008};
Plane Surface(1012) = {1011};
epaisseur=0.04;
//out[] = Extrude{0,0,epaisseur}{ Surface{1012}; };
Extrude{0,0,epaisseur}{ Surface{1012}; }
Transfinite Line {1007,1008} = ntheta;
Transfinite Line {1009,1010} = nr;
Transfinite Line {1014,1016} = nr;
Transfinite Line {1015,1017} = ntheta;
Transfinite Line {1019,1020} = nz;
Transfinite Line {1024,1028} = nz;
Transfinite Surface {1012} = {1002, 1003, 1004, 1005};
Transfinite Surface {1034} = {1006, 1007, 1012, 1016};
Transfinite Surface {1033} = {1002, 1005, 1006, 1016};
Transfinite Surface {1029} = {1004, 1005, 1012, 1016};
Transfinite Surface {1025} = {1003, 1004, 1007, 1012};
Transfinite Surface {1021} = {1002, 1003, 1006, 1007};
Recombine Surface{1012,1034,1033,1029,1025,1021}; // pourobtenir des quadrangles
Physical Surface(201) = {1021}; // bottom
Physical Surface(101) = {1029}; // top
Physical Surface(240) = {1025}; // r ext
Physical Surface(241) = {1033}; // r int
Physical Surface(52) = {1034,1012}; // interfaces
Physical Volume(1) = {1};
Transfinite Volume {1} = {1002, 1003, 1004, 1005, 1006, 1007, 1012, 1016};
# $Id: Makefile 16563 2008-08-15 16:08:28Z bangerth $
#target = EQUATIONS
target = laplacian
# For the small projects Makefile, you basically need to fill in only
# four fields.
#
# The first is the name of the application. It is assumed that the
# application name is the same as the base file name of the single C++
# file from which the application is generated.
# The second field determines whether you want to run your program in
# debug or optimized mode. The latter is significantly faster, but no
# run-time checking of parameters and internal states is performed, so
# you should set this value to `on' while you develop your program,
# and to `off' when running production computations.
debug-mode = off
# As third field, we need to give the path to the top-level deal.II
# directory. You need to adjust this to your needs. Since this path is
# probably the most often needed one in the Makefile internals, it is
# designated by a single-character variable, since that can be
# reference using /usr/share/deal.II only, i.e. without the parentheses that are
# required for most other parameters, as e.g. in $(target).
# The last field specifies the names of data and other files that
# shall be deleted when calling `make clean'. Object and backup files,
# executables and the like are removed anyway. Here, we give a list of
# files in the various output formats that deal.II supports.
clean-up-files = *gmv *gnuplot *gpl *eps *pov
#
#
# Usually, you will not need to change anything beyond this point.
#
#
# The next statement tell the `make' program where to find the
# deal.II top level directory and to include the file with the global
# settings
include /usr/share/deal.II/Make.global_options
# Since the whole project consists of only one file, we need not
# consider difficult dependencies. We only have to declare the
# libraries which we want to link to the object file, and there need
# to be two sets of libraries: one for the debug mode version of the
# application and one for the optimized mode. Here we have selected
# the versions for 2d. Note that the order in which the libraries are
# given here is important and that your applications won't link
# properly if they are given in another order.
#
# You may need to augment the lists of libraries when compiling your
# program for other dimensions, or when using third party libraries
libs.g = $(lib-deal2-2d.g) \
$(lib-deal2-3d.g) \
$(lib-lac.g) \
$(lib-base.g)
libs.o = $(lib-deal2-2d.o) \
$(lib-deal2-3d.o) \
$(lib-lac.o) \
$(lib-base.o)
# We now use the variable defined above which switch between debug and
# optimized mode to select the set of libraries to link with. Included
# in the list of libraries is the name of the object file which we
# will produce from the single C++ file. Note that by default we use
# the extension .g.o for object files compiled in debug mode and .o for
# object files in optimized mode (or whatever the local default on your
# system is instead of .o).
ifeq ($(debug-mode),on)
libraries = $(target).g.$(OBJEXT) $(libs.g)
else
libraries = $(target).$(OBJEXT) $(libs.o)
endif
# Now comes the first production rule: how to link the single object
# file produced from the single C++ file into the executable. Since
# this is the first rule in the Makefile, it is the one `make' selects
# if you call it without arguments.
$(target) : $(libraries)
@echo ============================ Linking $@
@$(CXX) -o $...@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
# To make running the application somewhat independent of the actual
# program name, we usually declare a rule `run' which simply runs the
# program. You can then run it by typing `make run'. This is also
# useful if you want to call the executable with arguments which do
# not change frequently. You may then want to add them to the
# following rule:
run: $(target)
@echo ============================ Running $<
@./$(target)$(EXEEXT)
# As a last rule to the `make' program, we define what to do when
# cleaning up a directory. This usually involves deleting object files
# and other automatically created files such as the executable itself,
# backup files, and data files. Since the latter are not usually quite
# diverse, you needed to declare them at the top of this file.
clean:
-rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
# Since we have not yet stated how to make an object file from a C++
# file, we should do so now. Since the many flags passed to the
# compiler are usually not of much interest, we suppress the actual
# command line using the `at' sign in the first column of the rules
# and write the string indicating what we do instead.
./%.g.$(OBJEXT) :
@echo ==============debug========= $(<F)
@$(CXX) $(CXXFLAGS.g) -c $< -o $@
./%.$(OBJEXT) :
@echo ==============optimized===== $(<F)
@$(CXX) $(CXXFLAGS.o) -c $< -o $@
# The following statement tells make that the rules `run' and `clean'
# are not expected to produce files of the same name as Makefile rules
# usually do.
.PHONY: run clean
# Finally there is a rule which you normally need not care much about:
# since the executable depends on some include files from the library,
# besides the C++ application file of course, it is necessary to
# re-generate the executable when one of the files it depends on has
# changed. The following rule to created a dependency file
# `Makefile.dep', which `make' uses to determine when to regenerate
# the executable. This file is automagically remade whenever needed,
# i.e. whenever one of the cc-/h-files changed. Make detects whether
# to remake this file upon inclusion at the bottom of this file.
#
# If the creation of Makefile.dep fails, blow it away and fail
Makefile.dep: $(target).cc Makefile \
$(shell echo /usr/include/deal.II/*/*.h)
@echo ============================ Remaking $@
@/usr/share/deal.II/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \
> $@ \
|| (rm -f $@ ; false)
@if test -s $@ ; then : else rm $@ ; fi
# To make the dependencies known to `make', we finally have to include
# them:
include Makefile.dep
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
