]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implementation of new laplace_transformation function.
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 18 May 2001 17:21:08 +0000 (17:21 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 18 May 2001 17:21:08 +0000 (17:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@4675 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_generator.h
deal.II/deal.II/source/grid/grid_generator.cc

index d64b3ddec14d8890ca6deb0465c5674e186d1cec..2ce974d7034d30ac44f67c0fa106fde4d423803a 100644 (file)
 
 
 #include <base/exceptions.h>
-
 template <int dim> class Point;
 template <int dim> class Triangulation;
+template <typename number> class Vector;
+template <typename number> class SparseMatrix;
+
 
+#include <map>
 
 /**
  * This class offers triangulations of some standard domains such as hypercubes,
@@ -95,7 +98,7 @@ template <int dim> 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 <int dim>
+    static void laplace_transformation (Triangulation<dim> &tria,
+                                       const std::map<unsigned int,Point<dim> > &new_points);
     
                                     /**
                                      * Exception
@@ -317,6 +344,17 @@ class GridGenerator
                                      */
     template <int dim>
     static void colorize_hyper_rectangle (Triangulation<dim> &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<double> &S,
+                              const std::map<unsigned int,double> &m,
+                              Vector<double> &u);
 };
 
 
index 78c2ed18b4a91c27439152f15f606712b55e7ba5..a3d30d1967e3f7f0f80860ba7352542a005fa868 100644 (file)
 //----------------------------  grid_generator.cc  ---------------------------
 
 
+#include <base/quadrature_lib.h>
+#include <base/thread_management.h>
+#include <lac/vector.h>
+#include <lac/vector_memory.h>
+#include <lac/filtered_matrix.h>
+#include <lac/precondition.h>
+#include <lac/solver_cg.h>
+#include <lac/sparse_matrix.h>
 #include <grid/grid_generator.h>
 #include <grid/tria.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/mapping_q1.h>
+#include <fe/fe_q.h>
+#include <numerics/matrices.h>
+
 #include <cmath>
 
 
@@ -723,6 +738,104 @@ void GridGenerator::hyper_shell (Triangulation<3>   &,
 #endif
 
 
+template <int dim>
+void GridGenerator::laplace_transformation (Triangulation<dim> &tria,
+                                           const std::map<unsigned int,Point<dim> > &new_points)
+{
+                                  // first provide everything that is
+                                  // needed for solving a Laplace
+                                  // equation.  
+  MappingQ1<dim> mapping_q1;  
+  FE_Q<dim> q1(1);
+  
+  DoFHandler<dim> 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<double> S(sparsity_pattern);
+  
+  QGauss4<dim> quadrature;
+  
+  MatrixCreator<dim>::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S);
+
+                                  // set up the boundary values for
+                                  // the laplace problem
+  std::vector<std::map<unsigned int,double> > m(dim);
+  typename std::map<unsigned int,Point<dim> >::const_iterator map_iter;
+
+                                  // fill these maps using the data
+                                  // given by new_points
+  typename DoFHandler<dim>::cell_iterator cell=dof_handler.begin_active(),
+                                endc=dof_handler.end();
+  typename DoFHandler<dim>::face_iterator face;
+  for (; cell!=endc; ++cell)
+    {
+      if (cell->at_boundary())
+       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+         {
+           face=cell->face(face_no);
+           if (face->at_boundary())
+             for (unsigned int vertex_no=0;
+                  vertex_no<GeometryInfo<dim>::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<dim; ++i)
+                   m[i].insert(pair<unsigned int,double> (
+                     face->vertex_dof_index(vertex_no, 0), map_iter->second(i)));
+               }
+         }
+    }
+  
+                                  // solve the dim problems with
+                                  // different right hand sides.
+  std::vector<Vector<double> > us(dim, Vector<double> (dof_handler.n_dofs()));
+  
+                                  // solve linear systems in parallel
+  Threads::ThreadManager thread_manager;
+  for (unsigned int i=0; i<dim; ++i)
+    Threads::spawn (thread_manager,
+                   Threads::encapsulate (&GridGenerator::laplace_solve)
+                   .collect_args (S, m[i], us[i]));
+  thread_manager.wait ();
+  
+                                  // change the coordinates of the
+                                  // points of the triangulation
+                                  // according to the computed values
+  for (cell=dof_handler.begin_active(); cell!=endc; ++cell)
+    for (unsigned int vertex_no=0;
+        vertex_no<GeometryInfo<dim>::vertices_per_cell; ++vertex_no)
+      {
+       Point<dim> &v=cell->vertex(vertex_no);
+       const unsigned int dof_index=cell->vertex_dof_index(vertex_no, 0);
+       for (unsigned int i=0; i<dim; ++i)
+         v(i)=us[i](dof_index);
+      }
+}
+
+
+void GridGenerator::laplace_solve (const SparseMatrix<double> &S,
+                                  const std::map<unsigned int,double> &m,
+                                  Vector<double> &u)
+{
+  const unsigned int n_dofs=S.n();
+  FilteredMatrix<SparseMatrix<double> > SF (S);
+  SolverControl control (1000, 1.e-10, false, false);
+  PrimitiveVectorMemory<Vector<double> > mem;
+  SolverCG<Vector<double> > solver (control, mem);
+  PreconditionJacobi<FilteredMatrix<SparseMatrix<double> > > prec;
+  Vector<double> 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<deal_II_dimension> &,
@@ -734,3 +847,7 @@ GridGenerator::hyper_cube (Triangulation<deal_II_dimension> &,
                           const double,
                           const double);
 
+template void
+GridGenerator::laplace_transformation (Triangulation<deal_II_dimension> &,
+                                      const std::map<unsigned int,Point<deal_II_dimension> > &);
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.