From 05b85ece7c5b8dfe438443d60553c998d27c9d33 Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Fri, 18 May 2001 17:21:08 +0000 Subject: [PATCH] Implementation of new laplace_transformation function. git-svn-id: https://svn.dealii.org/trunk@4675 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_generator.h | 42 ++++++- deal.II/deal.II/source/grid/grid_generator.cc | 117 ++++++++++++++++++ 2 files changed, 157 insertions(+), 2 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index d64b3ddec1..2ce974d703 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -15,10 +15,13 @@ #include - template class Point; template class Triangulation; +template class Vector; +template class SparseMatrix; + +#include /** * This class offers triangulations of some standard domains such as hypercubes, @@ -95,7 +98,7 @@ template class Triangulation; * of the eight child-cubes from one of their neighbors. * @end{itemize} * - * @author Wolfgang Bangerth, 1998, 1999. Slit domain by Stefan Nauber, 1999 + * @author Wolfgang Bangerth, 1998, 1999. Slit domain by Stefan Nauber, 1999. Laplace transformation by Ralf Hartmann, 2001. */ class GridGenerator { @@ -302,6 +305,30 @@ class GridGenerator const double inner_radius, const double outer_radius, const unsigned int n_cells = 0); + + + /** + * This function transformes the + * @p{Triangulation} @p{tria} + * smoothly to a domain that is + * described by the boundary + * points in the map + * @p{new_points}. This map maps + * the point indices to the + * boundary points in the + * transformed domain. + * + * Note, that the + * @p{Triangulation} is changed + * in-place, therefore you don't + * need to keep two + * triangulations, but the given + * triangulation is changed + * (overwritten). + */ + template + static void laplace_transformation (Triangulation &tria, + const std::map > &new_points); /** * Exception @@ -317,6 +344,17 @@ class GridGenerator */ template static void colorize_hyper_rectangle (Triangulation &tria); + + /** + * Multiply called by the + * @p{laplace_transformation} + * function. Solves the Laplace + * equation for one of the + * @p{dim} dimension. + */ + static void laplace_solve (const SparseMatrix &S, + const std::map &m, + Vector &u); }; diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index 78c2ed18b4..a3d30d1967 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -12,10 +12,25 @@ //---------------------------- grid_generator.cc --------------------------- +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include #include +#include +#include +#include +#include +#include +#include + #include @@ -723,6 +738,104 @@ void GridGenerator::hyper_shell (Triangulation<3> &, #endif +template +void GridGenerator::laplace_transformation (Triangulation &tria, + const std::map > &new_points) +{ + // first provide everything that is + // needed for solving a Laplace + // equation. + MappingQ1 mapping_q1; + FE_Q q1(1); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(q1); + SparsityPattern sparsity_pattern (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 (); + + SparseMatrix S(sparsity_pattern); + + QGauss4 quadrature; + + MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S); + + // set up the boundary values for + // the laplace problem + std::vector > m(dim); + typename std::map >::const_iterator map_iter; + + // fill these maps using the data + // given by new_points + typename DoFHandler::cell_iterator cell=dof_handler.begin_active(), + endc=dof_handler.end(); + typename DoFHandler::face_iterator face; + for (; cell!=endc; ++cell) + { + if (cell->at_boundary()) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face=cell->face(face_no); + if (face->at_boundary()) + for (unsigned int vertex_no=0; + vertex_no::vertices_per_face; ++vertex_no) + { + const unsigned int vertex_index=face->vertex_index(vertex_no); + map_iter=new_points.find(vertex_index); + for (unsigned int i=0; i ( + face->vertex_dof_index(vertex_no, 0), map_iter->second(i))); + } + } + } + + // solve the dim problems with + // different right hand sides. + std::vector > us(dim, Vector (dof_handler.n_dofs())); + + // solve linear systems in parallel + Threads::ThreadManager thread_manager; + for (unsigned int i=0; i::vertices_per_cell; ++vertex_no) + { + Point &v=cell->vertex(vertex_no); + const unsigned int dof_index=cell->vertex_dof_index(vertex_no, 0); + for (unsigned int i=0; i &S, + const std::map &m, + Vector &u) +{ + const unsigned int n_dofs=S.n(); + FilteredMatrix > SF (S); + SolverControl control (1000, 1.e-10, false, false); + PrimitiveVectorMemory > mem; + SolverCG > solver (control, mem); + PreconditionJacobi > > prec; + Vector f(n_dofs); + + SF.add_constraints(m); + prec.initialize (SF); + SF.apply_constraints (f, true); + solver.solve(SF, u, f, prec); +} + + // explicit instantiations template void GridGenerator::hyper_rectangle (Triangulation &, @@ -734,3 +847,7 @@ GridGenerator::hyper_cube (Triangulation &, const double, const double); +template void +GridGenerator::laplace_transformation (Triangulation &, + const std::map > &); + -- 2.39.5