From 05b85ece7c5b8dfe438443d60553c998d27c9d33 Mon Sep 17 00:00:00 2001
From: Ralf Hartmann <Ralf.Hartmann@dlr.de>
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 <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);
 };
 
 
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 <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> > &);
+
-- 
2.39.5