]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make it run. Not tested yet.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Aug 2011 18:58:05 +0000 (18:58 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Aug 2011 18:58:05 +0000 (18:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@24022 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-47/step-47.cc

index fadbfeed4c33b904480958f1734b7b817edcf177..420b63728ab72dc7c65b68e1c87b53374ae23b17 100644 (file)
@@ -53,6 +53,17 @@ using namespace dealii;
 
 
 
+double sign (double d)
+{
+  if (d > 0)
+    return 1;
+  else if (d < 0)
+    return -1;
+  else
+    return 0;
+}
+
+
 template <int dim>
 class LaplaceProblem
 {
@@ -63,12 +74,15 @@ class LaplaceProblem
     void run ();
 
   private:
+    double level_set (const Point<dim> &p) const;
+    Tensor<1,dim> grad_level_set (const Point<dim> &p) const;
+
     bool interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const;
-               std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
-               void append_quadrature(const Quadrature<dim> &plain_quadrature, 
-                                                  const std::vector<Point<dim> > &v      , 
-                                                                                                        std::vector<Point<dim> > &xfem_points, 
-                                                                                                        std::vector<double>      &xfem_weights);
+    std::pair<unsigned int, Quadrature<dim> > compute_quadrature(const Quadrature<dim> &plain_quadrature, const typename hp::DoFHandler<dim>::active_cell_iterator &cell, const std::vector<double> &level_set_values);
+    void append_quadrature(const Quadrature<dim> &plain_quadrature,
+                          const std::vector<Point<dim> > &v      ,
+                          std::vector<Point<dim> > &xfem_points,
+                          std::vector<double>      &xfem_weights);
 
     void setup_system ();
     void assemble_system ();
@@ -167,11 +181,38 @@ LaplaceProblem<dim>::~LaplaceProblem ()
 
 
 
+template <int dim>
+double
+LaplaceProblem<dim>::
+level_set (const Point<dim> &p) const
+{
+  return p.norm() - 0.5;
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+LaplaceProblem<dim>::
+grad_level_set (const Point<dim> &p) const
+{
+  return p / p.norm();
+}
+
+
+
 template <int dim>
 bool
 LaplaceProblem<dim>::
 interface_intersects_cell (const typename Triangulation<dim>::cell_iterator &cell) const
 {
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell-1; ++v)
+    if (level_set(cell->vertex(v)) * level_set(cell->vertex(v+1)) < 0)
+      return true;
+
+                                  // we get here only if all vertices
+                                  // have the same sign, which means
+                                  // that the cell is not intersected
   return false;
 }
 
@@ -230,9 +271,9 @@ void LaplaceProblem<dim>::setup_system ()
 
 
   constraints.clear ();
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          constraints);
-  constraints.close ();
+//TODO: fix this, it currently crashes
+  // DoFTools::make_hanging_node_constraints (dof_handler,
+  //                                      constraints);
 
                                   // now constrain those enriched
                                   // DoFs that are on cells that are
@@ -252,7 +293,9 @@ void LaplaceProblem<dim>::setup_system ()
                                       // constrain these DoFs
       for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
        if (vertex_is_on_intersected_cell[cell->vertex_index(v)] == false)
-         constraints.add_line (cell->vertex_dof_index(v,1));
+         constraints.add_line (cell->vertex_dof_index(v,1,cell->active_fe_index()));
+
+//TODO: component 1 must satisfy zero boundary conditions
   constraints.close();
 
 
@@ -275,6 +318,9 @@ void LaplaceProblem<dim>::assemble_system ()
   FEValues<dim> plain_fe_values (fe_collection[0], quadrature_formula,
                                 update_values    |  update_gradients |
                                 update_quadrature_points  |  update_JxW_values);
+  FEValues<dim> enriched_fe_values (fe_collection[1], quadrature_formula,
+                                            update_values    |  update_gradients |
+                                            update_quadrature_points  |  update_JxW_values);
 
   const unsigned int   n_q_points    = quadrature_formula.size();
 
@@ -290,30 +336,13 @@ void LaplaceProblem<dim>::assemble_system ()
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
 
-       std::vector<double> level_set_values;
-       level_set_values.push_back(1.);
-       level_set_values.push_back(1.);
-       level_set_values.push_back(1.);
-       level_set_values.push_back(-3.);
+  std::vector<double> level_set_values;
+  level_set_values.push_back(1.);
+  level_set_values.push_back(1.);
+  level_set_values.push_back(1.);
+  level_set_values.push_back(-3.);
   for (; cell!=endc; ++cell)
     {
-                       std::pair<unsigned int, Quadrature<dim> > type_and_quadrature = compute_quadrature(quadrature_formula, cell, level_set_values);
-
-                       std::cout << "type : " << type_and_quadrature.first << std::endl;
-                       std::vector<Point<dim> > points = type_and_quadrature.second.get_points();
-                       std::vector<double> weights = type_and_quadrature.second.get_weights();
-                       std::string filename = "points.dat";
-                       std::ofstream output (filename.c_str());
-                       output << "#xfem quadrature Points" << std::endl;
-                       for (unsigned int i=0; i<points.size(); i++)
-                               output << points[i] << std::endl;
-                       std::string filename2 = "weights.dat";
-                       std::ofstream output2 (filename2.c_str());
-                       output2 << "#xfem Weights" << std::endl;
-                       for (unsigned int i=0; i<weights.size(); i++)
-                               output2 << weights[i] << std::endl;
-                       assert(0);
-
       const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
       cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
       cell_rhs.reinit (dofs_per_cell);
@@ -321,24 +350,121 @@ void LaplaceProblem<dim>::assemble_system ()
       cell_matrix = 0;
       cell_rhs = 0;
 
-      plain_fe_values.reinit (cell);
+      if (cell->active_fe_index() == 0)
+       {
+         plain_fe_values.reinit (cell);
 
-      coefficient.value_list (plain_fe_values.get_quadrature_points(),
-                             coefficient_values);
+         coefficient_values.resize (plain_fe_values.n_quadrature_points);
+         coefficient.value_list (plain_fe_values.get_quadrature_points(),
+                                 coefficient_values);
 
-      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-       for (unsigned int i=0; i<dofs_per_cell; ++i)
-         {
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += (coefficient_values[q_point] *
-                                  plain_fe_values.shape_grad(i,q_point) *
-                                  plain_fe_values.shape_grad(j,q_point) *
-                                  plain_fe_values.JxW(q_point));
-
-           cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
-                           1.0 *
-                           plain_fe_values.JxW(q_point));
-         }
+         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             {
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
+                 cell_matrix(i,j) += (coefficient_values[q_point] *
+                                      plain_fe_values.shape_grad(i,q_point) *
+                                      plain_fe_values.shape_grad(j,q_point) *
+                                      plain_fe_values.JxW(q_point));
+
+               cell_rhs(i) += (plain_fe_values.shape_value(i,q_point) *
+                               1.0 *
+                               plain_fe_values.JxW(q_point));
+             }
+       }
+      else
+       {
+         Assert (cell->active_fe_index() == 1, ExcInternalError());
+
+         std::auto_ptr<FEValues<dim> > custom_fe_values;
+
+         if (interface_intersects_cell(cell) == false)
+           enriched_fe_values.reinit (cell);
+         else
+           {
+             std::vector<double> level_set_values (GeometryInfo<dim>::vertices_per_cell);
+             for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+               level_set_values[v] = level_set (cell->vertex(v));
+
+             custom_fe_values
+               .reset(new FEValues<dim> (fe_collection[1],
+                                         compute_quadrature(quadrature_formula, cell,
+                                                            level_set_values).second,
+                                         update_values    |  update_gradients |
+                                         update_quadrature_points  |  update_JxW_values));
+             custom_fe_values->reinit (cell);
+           }
+
+         FEValues<dim> &this_fe_values = ((interface_intersects_cell(cell) == false)
+                                          ?
+                                          enriched_fe_values
+                                          :
+                                          *custom_fe_values);
+
+         coefficient_values.resize (this_fe_values.n_quadrature_points);
+         coefficient.value_list (this_fe_values.get_quadrature_points(),
+                                 coefficient_values);
+
+         for (unsigned int q_point=0; q_point<this_fe_values.n_quadrature_points; ++q_point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             if (cell->get_fe().system_to_component_index(i).first == 0)
+               {
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   if (cell->get_fe().system_to_component_index(j).first == 0)
+                     cell_matrix(i,j) += (coefficient_values[q_point] *
+                                          this_fe_values.shape_grad(i,q_point) *
+                                          this_fe_values.shape_grad(j,q_point) *
+                                          this_fe_values.JxW(q_point));
+                   else
+                     cell_matrix(i,j) += (coefficient_values[q_point] *
+                                          this_fe_values.shape_grad(i,q_point) *
+                                          (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_grad(j,q_point)
+                                           +
+                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_value(j,q_point)) *
+                                          this_fe_values.JxW(q_point));
+
+                 cell_rhs(i) += (this_fe_values.shape_value(i,q_point) *
+                                 1.0 *
+                                 this_fe_values.JxW(q_point));
+               }
+             else
+               {
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   if (cell->get_fe().system_to_component_index(j).first == 0)
+                     cell_matrix(i,j) += (coefficient_values[q_point] *
+                                          (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_grad(i,q_point)
+                                           +
+                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_value(i,q_point)) *
+                                          this_fe_values.shape_grad(j,q_point) *
+                                          this_fe_values.JxW(q_point));
+                   else
+                     cell_matrix(i,j) += (coefficient_values[q_point] *
+                                          (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_grad(i,q_point)
+                                           +
+                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_value(i,q_point)) *
+                                          (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_grad(j,q_point)
+                                           +
+                                           grad_level_set(this_fe_values.quadrature_point(q_point)) *
+                                           sign(level_set(this_fe_values.quadrature_point(q_point))) *
+                                           this_fe_values.shape_value(j,q_point)) *
+                                          this_fe_values.JxW(q_point));
+
+                 cell_rhs(i) += (std::fabs(level_set(this_fe_values.quadrature_point(q_point))) *
+                                 this_fe_values.shape_value(i,q_point) *
+                                 1.0 *
+                                 this_fe_values.JxW(q_point));
+               }
+       }
 
       local_dof_indices.resize (dofs_per_cell);
       cell->get_dof_indices (local_dof_indices);
@@ -368,408 +494,410 @@ void LaplaceProblem<dim>::assemble_system ()
 // are considered.
 // Type 1: there is not cut. Type 2: a corner of the element is cut. Type 3: two corners are cut.
 
-       template <int dim>
-       std::pair<unsigned int, Quadrature<dim> >
-       LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
-               const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
-               const std::vector<double> &level_set_values                    )
+template <int dim>
+std::pair<unsigned int, Quadrature<dim> >
+LaplaceProblem<dim>::compute_quadrature (const Quadrature<dim> &plain_quadrature,
+                                        const typename hp::DoFHandler<dim>::active_cell_iterator &cell,
+                                        const std::vector<double> &level_set_values                    )
 {
 
-       unsigned int type = 0;
+  unsigned int type = 0;
 
-       // find the type of cut
-       int sign_ls[GeometryInfo<dim>::vertices_per_cell];
-       for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-       {
-               if (level_set_values[v] > 0) sign_ls[v] = 1;
-               else if (level_set_values[v] < 0) sign_ls[v] = -1;
-               else sign_ls[v] = 0;
-       }
+                                  // find the type of cut
+  int sign_ls[GeometryInfo<dim>::vertices_per_cell];
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+    {
+      if (level_set_values[v] > 0) sign_ls[v] = 1;
+      else if (level_set_values[v] < 0) sign_ls[v] = -1;
+      else sign_ls[v] = 0;
+    }
 
-       // the sign of the level set function at the 4 nodes of the elements can be positive + or negative -
-       // depending on the sign of the level set function we have the folloing three classes of decomposition
-  // type 1: ++++, ----
-       // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+
-       // type 3: +--+, ++--, +-+-, -++-, --++, -+-+
+                                  // the sign of the level set function at the 4 nodes of the elements can be positive + or negative -
+                                  // depending on the sign of the level set function we have the folloing three classes of decomposition
+                                  // type 1: ++++, ----
+                                  // type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+
+                                  // type 3: +--+, ++--, +-+-, -++-, --++, -+-+
 
-       if ( sign_ls[0]==sign_ls[1] & sign_ls[0]==sign_ls[2] & sign_ls[0]==sign_ls[3] ) type =1;
-       else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 ) type = 2;
-       else type = 3;
+  if ( sign_ls[0]==sign_ls[1] & sign_ls[0]==sign_ls[2] & sign_ls[0]==sign_ls[3] ) type =1;
+  else if ( sign_ls[0]*sign_ls[1]*sign_ls[2]*sign_ls[3] < 0 ) type = 2;
+  else type = 3;
 
-       unsigned int Pos = 100;
+  unsigned int Pos = 100;
 
-       Point<dim> v0(0,0);
-       Point<dim> v1(1,0);
-       Point<dim> v2(0,1);
-       Point<dim> v3(1,1);
+  Point<dim> v0(0,0);
+  Point<dim> v1(1,0);
+  Point<dim> v2(0,1);
+  Point<dim> v3(1,1);
 
-       Point<dim> A(0,0);
-       Point<dim> B(0,0);
-       Point<dim> C(0,0);
-       Point<dim> D(0,0);
-       Point<dim> E(0,0);
-       Point<dim> F(0,0);
+  Point<dim> A(0,0);
+  Point<dim> B(0,0);
+  Point<dim> C(0,0);
+  Point<dim> D(0,0);
+  Point<dim> E(0,0);
+  Point<dim> F(0,0);
 
-       if (type == 1)
-         return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
+  if (type == 1)
+    return std::pair<unsigned int, Quadrature<dim> >(1, plain_quadrature);
 
-       if (type==2)
-       {
-               const unsigned int   n_q_points    = plain_quadrature.size();
+  if (type==2)
+    {
+      const unsigned int   n_q_points    = plain_quadrature.size();
 
-               // loop over all subelements for integration
-               // in type 2 there are 5 subelements
+                                      // loop over all subelements for integration
+                                      // in type 2 there are 5 subelements
 
-               Quadrature<dim> xfem_quadrature(5*n_q_points);
+      Quadrature<dim> xfem_quadrature(5*n_q_points);
 
-               std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
+      std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
 
-               if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0;
-    else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1;
-    else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2;
-    else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3;
-               else assert(0); // error message
+      if (sign_ls[0]!=sign_ls[1] && sign_ls[0]!=sign_ls[2] && sign_ls[0]!=sign_ls[3]) Pos = 0;
+      else if (sign_ls[1]!=sign_ls[0] && sign_ls[1]!=sign_ls[2] && sign_ls[1]!=sign_ls[3]) Pos = 1;
+      else if (sign_ls[2]!=sign_ls[0] && sign_ls[2]!=sign_ls[1] && sign_ls[2]!=sign_ls[3]) Pos = 2;
+      else if (sign_ls[3]!=sign_ls[0] && sign_ls[3]!=sign_ls[1] && sign_ls[3]!=sign_ls[2]) Pos = 3;
+      else assert(0); // error message
 
-               // Find cut coordinates
+                                      // Find cut coordinates
 
-               // deal.ii local coordinates
+                                      // deal.ii local coordinates
 
-               //    2-------3
-               //    |       |
-               //              |       |
-               //              |       |
-               //              0-------1
+                                      //    2-------3
+                                      //    |       |
+                                      //               |       |
+                                      //               |       |
+                                      //               0-------1
 
-               if (Pos == 0)
-               {
-                       A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]);
-                       B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]);
-                       A(1) = 0.;
-                       B(0) = 0.;
-                       C(0) = 0.5*( A(0) + B(0) );
-                       C(1) = 0.5*( A(1) + B(1) );
-                       D(0) = 2./3. * C(0);
-                       D(1) = 2./3. * C(1);
-                       E(0) = 0.5*A(0);
-                       E(1) = 0.;
-                       F(0) = 0.;
-                       F(1) = 0.5*B(1);
-               }
-               else if (Pos == 1)
-               {
-                       A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]);
-                       B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]);
-                       A(1) = 0.;
-                       B(0) = 1.;
-                       C(0) = 0.5*( A(0) + B(0) );
-                       C(1) = 0.5*( A(1) + B(1) );
-                       D(0) = 1./3. + 2./3. * C(0);
-                       D(1) = 2./3. * C(1);
-                       E(0) = 0.5*(1 + A(0));
-                       E(1) = 0.;
-                       F(0) = 1.;
-                       F(1) = 0.5*B(1);
-                       }
-               else if (Pos == 2)
-               {
-                       A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]);
-                       B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]);
-                       A(1) = 1.;
-                       B(0) = 0.;
-                       C(0) = 0.5*( A(0) + B(0) );
-                       C(1) = 0.5*( A(1) + B(1) );
-                       D(0) = 2./3. * C(0);
-                       D(1) = 1./3. + 2./3. * C(1);
-                       E(0) = 0.5* A(0);
-                       E(1) = 1.;
-                       F(0) = 0.;
-                       F(1) = 0.5*( 1. + B(1) );
-                       }
-               else if (Pos == 3)
-               {
-                       A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]);
-                       B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]);
-                       A(1) = 1.;
-                       B(0) = 1.;
-                       C(0) = 0.5*( A(0) + B(0) );
-                       C(1) = 0.5*( A(1) + B(1) );
-                       D(0) = 1./3. + 2./3. * C(0);
-                       D(1) = 1./3. + 2./3. * C(1);
-                       E(0) = 0.5*( 1. + A(0) );
-                       E(1) = 1.;
-                       F(0) = 1.;
-                       F(1) = 0.5*( 1. + B(1) );
-                       }
-
-               //std::cout << A << std::endl;
-               //std::cout << B << std::endl;
-               //std::cout << C << std::endl;
-               //std::cout << D << std::endl;
-               //std::cout << E << std::endl;
-               //std::cout << F << std::endl;
-
-               std::string filename = "vertices.dat";
-               std::ofstream output (filename.c_str());
-               output << "#vertices of xfem subcells" << std::endl;
-               output << v0(0) << "   " << v0(1) << std::endl;
-               output << v1(0) << "   " << v1(1) << std::endl;
-               output << v3(0) << "   " << v3(1) << std::endl;
-               output << v2(0) << "   " << v2(1) << std::endl;
-               output << std::endl;
-               output << A(0) << "   " << A(1) << std::endl;
-               output << B(0) << "   " << B(1) << std::endl;
-               output << std::endl;
-               output << C(0) << "   " << C(1) << std::endl;
-               output << D(0) << "   " << D(1) << std::endl;
-               output << std::endl;
-               output << D(0) << "   " << D(1) << std::endl;
-               output << E(0) << "   " << E(1) << std::endl;
-               output << std::endl;
-               output << D(0) << "   " << D(1) << std::endl;
-               output << F(0) << "   " << F(1) << std::endl;
-               output << std::endl;
-
-               if (Pos==0)
-                       output << v3(0) << "   " << v3(1) << std::endl;
-               else if (Pos==1)
-                       output << v2(0) << "   " << v2(1) << std::endl;
-               else if (Pos==2)
-                       output << v1(0) << "   " << v1(1) << std::endl;
-               else if (Pos==3)
-                       output << v0(0) << "   " << v0(1) << std::endl;
-               output << C(0) << "   " << C(1) << std::endl;
-
-               Point<dim> subcell_vertices[10];
-               subcell_vertices[0] = v0;
-               subcell_vertices[1] = v1;
-               subcell_vertices[2] = v2;
-               subcell_vertices[3] = v3;
-               subcell_vertices[4] = A;
-               subcell_vertices[5] = B;
-               subcell_vertices[6] = C;
-               subcell_vertices[7] = D;
-               subcell_vertices[8] = E;
-               subcell_vertices[9] = F;
-
-               std::vector<Point<dim> > xfem_points;
-               std::vector<double>      xfem_weights;
-
-               // lookup table for the decomposition
-
-               if (dim==2)
+      if (Pos == 0)
+       {
+         A[0] = 1. - level_set_values[1]/(level_set_values[1]-level_set_values[0]);
+         B[1] = 1. - level_set_values[2]/(level_set_values[2]-level_set_values[0]);
+         A(1) = 0.;
+         B(0) = 0.;
+         C(0) = 0.5*( A(0) + B(0) );
+         C(1) = 0.5*( A(1) + B(1) );
+         D(0) = 2./3. * C(0);
+         D(1) = 2./3. * C(1);
+         E(0) = 0.5*A(0);
+         E(1) = 0.;
+         F(0) = 0.;
+         F(1) = 0.5*B(1);
+       }
+      else if (Pos == 1)
+       {
+         A[0] = level_set_values[0]/(level_set_values[0]-level_set_values[1]);
+         B[1] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[1]);
+         A(1) = 0.;
+         B(0) = 1.;
+         C(0) = 0.5*( A(0) + B(0) );
+         C(1) = 0.5*( A(1) + B(1) );
+         D(0) = 1./3. + 2./3. * C(0);
+         D(1) = 2./3. * C(1);
+         E(0) = 0.5*(1 + A(0));
+         E(1) = 0.;
+         F(0) = 1.;
+         F(1) = 0.5*B(1);
+       }
+      else if (Pos == 2)
+       {
+         A[0] = 1 - level_set_values[3]/(level_set_values[3]-level_set_values[2]);
+         B[1] = level_set_values[0]/(level_set_values[0]-level_set_values[2]);
+         A(1) = 1.;
+         B(0) = 0.;
+         C(0) = 0.5*( A(0) + B(0) );
+         C(1) = 0.5*( A(1) + B(1) );
+         D(0) = 2./3. * C(0);
+         D(1) = 1./3. + 2./3. * C(1);
+         E(0) = 0.5* A(0);
+         E(1) = 1.;
+         F(0) = 0.;
+         F(1) = 0.5*( 1. + B(1) );
+       }
+      else if (Pos == 3)
+       {
+         A[0] = level_set_values[2]/(level_set_values[2]-level_set_values[3]);
+         B[1] = level_set_values[1]/(level_set_values[1]-level_set_values[3]);
+         A(1) = 1.;
+         B(0) = 1.;
+         C(0) = 0.5*( A(0) + B(0) );
+         C(1) = 0.5*( A(1) + B(1) );
+         D(0) = 1./3. + 2./3. * C(0);
+         D(1) = 1./3. + 2./3. * C(1);
+         E(0) = 0.5*( 1. + A(0) );
+         E(1) = 1.;
+         F(0) = 1.;
+         F(1) = 0.5*( 1. + B(1) );
+       }
+
+                                      //std::cout << A << std::endl;
+                                      //std::cout << B << std::endl;
+                                      //std::cout << C << std::endl;
+                                      //std::cout << D << std::endl;
+                                      //std::cout << E << std::endl;
+                                      //std::cout << F << std::endl;
+
+      std::string filename = "vertices.dat";
+      std::ofstream output (filename.c_str());
+      output << "#vertices of xfem subcells" << std::endl;
+      output << v0(0) << "   " << v0(1) << std::endl;
+      output << v1(0) << "   " << v1(1) << std::endl;
+      output << v3(0) << "   " << v3(1) << std::endl;
+      output << v2(0) << "   " << v2(1) << std::endl;
+      output << std::endl;
+      output << A(0) << "   " << A(1) << std::endl;
+      output << B(0) << "   " << B(1) << std::endl;
+      output << std::endl;
+      output << C(0) << "   " << C(1) << std::endl;
+      output << D(0) << "   " << D(1) << std::endl;
+      output << std::endl;
+      output << D(0) << "   " << D(1) << std::endl;
+      output << E(0) << "   " << E(1) << std::endl;
+      output << std::endl;
+      output << D(0) << "   " << D(1) << std::endl;
+      output << F(0) << "   " << F(1) << std::endl;
+      output << std::endl;
+
+      if (Pos==0)
+       output << v3(0) << "   " << v3(1) << std::endl;
+      else if (Pos==1)
+       output << v2(0) << "   " << v2(1) << std::endl;
+      else if (Pos==2)
+       output << v1(0) << "   " << v1(1) << std::endl;
+      else if (Pos==3)
+       output << v0(0) << "   " << v0(1) << std::endl;
+      output << C(0) << "   " << C(1) << std::endl;
+
+      Point<dim> subcell_vertices[10];
+      subcell_vertices[0] = v0;
+      subcell_vertices[1] = v1;
+      subcell_vertices[2] = v2;
+      subcell_vertices[3] = v3;
+      subcell_vertices[4] = A;
+      subcell_vertices[5] = B;
+      subcell_vertices[6] = C;
+      subcell_vertices[7] = D;
+      subcell_vertices[8] = E;
+      subcell_vertices[9] = F;
+
+      std::vector<Point<dim> > xfem_points;
+      std::vector<double>      xfem_weights;
+
+                                      // lookup table for the decomposition
+
+      if (dim==2)
+       {
+         unsigned int subcell_v_indices[4][5][4] = {
+               {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}},
+               {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}},
+               {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}},
+               {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}}
+         };
+
+         std::cout << "Pos       : " << Pos << std::endl;
+         for (unsigned int subcell = 0; subcell<5; subcell++)
+           {
+                                              //std::cout << "subcell   : " << subcell << std::endl;
+             std::vector<Point<dim> > vertices;
+             for (unsigned int i=0; i<4; i++)
                {
-                        unsigned int subcell_v_indices[4][5][4] = {
-                                {{0,8,9,7}, {9,7,5,6}, {8,4,7,6}, {5,6,2,3}, {6,4,3,1}},
-                                {{8,1,7,9}, {4,8,6,7}, {6,7,5,9}, {0,4,2,6}, {2,6,3,5}},
-                                {{9,7,2,8}, {5,6,9,7}, {6,4,7,8}, {0,1,5,6}, {6,1,4,3}},
-                                {{7,9,8,3}, {4,6,8,7}, {6,5,7,9}, {0,6,2,4}, {0,1,6,5}}
-                        };
-
-                        std::cout << "Pos       : " << Pos << std::endl;
-                        for (unsigned int subcell = 0; subcell<5; subcell++)
-                        {
-                                //std::cout << "subcell   : " << subcell << std::endl;
-                                std::vector<Point<dim> > vertices;
-                                for (unsigned int i=0; i<4; i++)
-                                {
-                                        vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
-                                        //std::cout << "i         : " << i << std::endl;
-                                        //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
-                                        //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
-                                }
-                                //std::cout << std::endl;
-                                // create quadrature rule
-                                append_quadrature( plain_quadrature,
-                                                                     vertices,
-                                                                                                               xfem_points,
-                                                                                                               xfem_weights);
-                                //initialize xfem_quadrature with quadrature points of all subelements
-                                xfem_quadrature.initialize(xfem_points, xfem_weights);
-                        }
+                 vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
+                                                  //std::cout << "i         : " << i << std::endl;
+                                                  //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
+                                                  //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
                }
-
-               return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
+                                              //std::cout << std::endl;
+                                              // create quadrature rule
+             append_quadrature( plain_quadrature,
+                                vertices,
+                                xfem_points,
+                                xfem_weights);
+                                              //initialize xfem_quadrature with quadrature points of all subelements
+             xfem_quadrature.initialize(xfem_points, xfem_weights);
+           }
        }
 
-       // Type three decomposition
-       // (+--+, ++--, +-+-, -++-, --++, -+-+)
+      Assert (xfem_quadrature.size() == 20, ExcInternalError());
+      return std::pair<unsigned int, Quadrature<dim> >(2, xfem_quadrature);
+    }
 
-       if (type==3)
-       {
-               const unsigned int   n_q_points    = plain_quadrature.size();
+                                  // Type three decomposition
+                                  // (+--+, ++--, +-+-, -++-, --++, -+-+)
 
-               // loop over all subelements for integration
-               // in type 2 there are 5 subelements
+  if (type==3)
+    {
+      const unsigned int   n_q_points    = plain_quadrature.size();
 
-               Quadrature<dim> xfem_quadrature(5*n_q_points);
+                                      // loop over all subelements for integration
+                                      // in type 2 there are 5 subelements
 
-               std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
+      Quadrature<dim> xfem_quadrature(5*n_q_points);
 
-               if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] )
-               {
-                       Pos = 0;
-                       A(0) = 0.;
-                       A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2]));
-                       B(0) = 1.;
-                       B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3]));
-               }
-               else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] )
-               {
-                       Pos = 1;
-                       A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1]));
-                       A(1) = 0.;
-                       B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3]));
-                       B(1) = 1.;
-               }
-               else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] )
-               {
-                       std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl;
-                       assert(0);
-               }
-               else
-               {
-                       std::cout << "Error: the level set function has not the right values" << std::endl;
-                       assert(0);
-               }
+      std::vector<Point<dim> > v(GeometryInfo<dim>::vertices_per_cell);
 
-               //std::cout << "Pos " << Pos << std::endl;
-               //std::cout << A << std::endl;
-               //std::cout << B << std::endl;
-               std::string filename = "vertices.dat";
-               std::ofstream output (filename.c_str());
-               output << "#vertices of xfem subcells" << std::endl;
-               output << A(0) << "   " << A(1) << std::endl;
-               output << B(0) << "   " << B(1) << std::endl;
-
-               //fill xfem_quadrature
-               Point<dim> subcell_vertices[6];
-               subcell_vertices[0] = v0;
-               subcell_vertices[1] = v1;
-               subcell_vertices[2] = v2;
-               subcell_vertices[3] = v3;
-               subcell_vertices[4] = A;
-               subcell_vertices[5] = B;
-       
-               std::vector<Point<dim> > xfem_points;
-               std::vector<double>      xfem_weights;
-
-               if (dim==2)
+      if ( sign_ls[0]==sign_ls[1] && sign_ls[2]==sign_ls[3] )
+       {
+         Pos = 0;
+         A(0) = 0.;
+         A(1) = level_set_values[0]/((level_set_values[0]-level_set_values[2]));
+         B(0) = 1.;
+         B(1) = level_set_values[1]/((level_set_values[1]-level_set_values[3]));
+       }
+      else if ( sign_ls[0]==sign_ls[2] && sign_ls[1]==sign_ls[3] )
+       {
+         Pos = 1;
+         A(0) = level_set_values[0]/((level_set_values[0]-level_set_values[1]));
+         A(1) = 0.;
+         B(0) = level_set_values[2]/((level_set_values[2]-level_set_values[3]));
+         B(1) = 1.;
+       }
+      else if ( sign_ls[0]==sign_ls[3] && sign_ls[1]==sign_ls[2] )
+       {
+         std::cout << "Error: the element has two cut lines and this is not allowed" << std::endl;
+         assert(0);
+       }
+      else
+       {
+         std::cout << "Error: the level set function has not the right values" << std::endl;
+         assert(0);
+       }
+
+                                      //std::cout << "Pos " << Pos << std::endl;
+                                      //std::cout << A << std::endl;
+                                      //std::cout << B << std::endl;
+      std::string filename = "vertices.dat";
+      std::ofstream output (filename.c_str());
+      output << "#vertices of xfem subcells" << std::endl;
+      output << A(0) << "   " << A(1) << std::endl;
+      output << B(0) << "   " << B(1) << std::endl;
+
+                                      //fill xfem_quadrature
+      Point<dim> subcell_vertices[6];
+      subcell_vertices[0] = v0;
+      subcell_vertices[1] = v1;
+      subcell_vertices[2] = v2;
+      subcell_vertices[3] = v3;
+      subcell_vertices[4] = A;
+      subcell_vertices[5] = B;
+
+      std::vector<Point<dim> > xfem_points;
+      std::vector<double>      xfem_weights;
+
+      if (dim==2)
+       {
+         unsigned int subcell_v_indices[2][2][4] = {
+               {{0,1,4,5}, {4,5,2,3}},
+               {{0,4,2,5}, {4,1,5,3}}
+         };
+
+                                          //std::cout << "Pos       : " << Pos << std::endl;
+         for (unsigned int subcell = 0; subcell<2; subcell++)
+           {
+                                              //std::cout << "subcell   : " << subcell << std::endl;
+             std::vector<Point<dim> > vertices;
+             for (unsigned int i=0; i<4; i++)
                {
-                        unsigned int subcell_v_indices[2][2][4] = {
-                                {{0,1,4,5}, {4,5,2,3}},
-                                {{0,4,2,5}, {4,1,5,3}}
-                        };
-
-                        //std::cout << "Pos       : " << Pos << std::endl;
-                        for (unsigned int subcell = 0; subcell<2; subcell++)
-                        {
-                                //std::cout << "subcell   : " << subcell << std::endl;
-                                std::vector<Point<dim> > vertices;
-                                for (unsigned int i=0; i<4; i++)
-                                {
-                                        vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
-                                        //std::cout << "i         : " << i << std::endl;
-                                        //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
-                                        //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
-                                }
-                                //std::cout << std::endl;
-                                // create quadrature rule
-                                append_quadrature( plain_quadrature,
-                                                                     vertices,
-                                                                                                               xfem_points,
-                                                                                                               xfem_weights);
-                                //initialize xfem_quadrature with quadrature points of all subelements
-                                xfem_quadrature.initialize(xfem_points, xfem_weights);
-                        }
+                 vertices.push_back( subcell_vertices[subcell_v_indices[Pos][subcell][i]] );
+                                                  //std::cout << "i         : " << i << std::endl;
+                                                  //std::cout << "subcell v : " << subcell_v_indices[Pos][subcell][i] << std::endl;
+                                                  //std::cout << vertices[i](0) << "  " << vertices[i](1) << std::endl;
                }
-               return std::pair<unsigned int, Quadrature<dim> >(3, xfem_quadrature);
+                                              //std::cout << std::endl;
+                                              // create quadrature rule
+             append_quadrature( plain_quadrature,
+                                vertices,
+                                xfem_points,
+                                xfem_weights);
+                                              //initialize xfem_quadrature with quadrature points of all subelements
+             xfem_quadrature.initialize(xfem_points, xfem_weights);
+           }
        }
+      Assert (xfem_quadrature.size() == 8, ExcInternalError());
+      return std::pair<unsigned int, Quadrature<dim> >(3, xfem_quadrature);
+    }
 
-       return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
+  return std::pair<unsigned int, Quadrature<dim> >(0, plain_quadrature);;
 
 }
 
-       template <int dim>
+template <int dim>
 void LaplaceProblem<dim>::append_quadrature ( const Quadrature<dim> &plain_quadrature,
                                              const std::vector<Point<dim> > &v,
                                              std::vector<Point<dim> > &xfem_points,
                                              std::vector<double>      &xfem_weights)
 
 {
-       // Project integration points into sub-elements.
-       // This maps quadrature points from a reference element to a subelement of a reference element.
-       // To implement the action of this map the coordinates of the subelements have been calculated (A(0)...F(0),A(1)...F(1))
-       // the coordinates of the quadrature points are given by the bi-linear map defined by the form functions
-       // $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the $\phi_j$ are the shape functions of the FEQ.
+                                  // Project integration points into sub-elements.
+                                  // This maps quadrature points from a reference element to a subelement of a reference element.
+                                  // To implement the action of this map the coordinates of the subelements have been calculated (A(0)...F(0),A(1)...F(1))
+                                  // the coordinates of the quadrature points are given by the bi-linear map defined by the form functions
+                                  // $x^\prime_i = \sum_j v^\prime \phi_j(x^hat_i)$, where the $\phi_j$ are the shape functions of the FEQ.
 
-       unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
+  unsigned int n_v = GeometryInfo<dim>::vertices_per_cell;
 
-       std::vector<Point<dim> > q_points = plain_quadrature.get_points();
-       std::vector<Point<dim> > q_transf(q_points.size());
-       std::vector<double> W = plain_quadrature.get_weights();
-       std::vector<double> phi(n_v);
-       std::vector<Tensor<1,dim> > grad_phi(n_v);
+  std::vector<Point<dim> > q_points = plain_quadrature.get_points();
+  std::vector<Point<dim> > q_transf(q_points.size());
+  std::vector<double> W = plain_quadrature.get_weights();
+  std::vector<double> phi(n_v);
+  std::vector<Tensor<1,dim> > grad_phi(n_v);
 
-       const unsigned int   n_q_points    = plain_quadrature.size();
+  const unsigned int   n_q_points    = plain_quadrature.size();
 
-       std::vector<double> JxW(n_q_points);
+  std::vector<double> JxW(n_q_points);
 
-       for ( unsigned int i = 0; i < n_q_points; i++)
+  for ( unsigned int i = 0; i < n_q_points; i++)
+    {
+      switch (dim)
        {
-         switch (dim)
-           {
-             case 2:
-             {
-               double xi  = q_points[i](0);
-               double eta = q_points[i](1);
-
-                                                // Define shape functions on reference element
-                                                // we consider a bi-linear mapping
-               phi[0] = (1. - xi) * (1. - eta);
-               phi[1] = xi * (1. - eta);
-               phi[2] = (1. - xi) * eta;
-               phi[3] = xi * eta;
-
-               grad_phi[0][0] = (-1. + eta);
-               grad_phi[1][0] = (1. - eta);
-               grad_phi[2][0] = -eta;
-               grad_phi[3][0] = eta;
-
-               grad_phi[0][1] = (-1. + xi);
-               grad_phi[1][1] = -xi;
-               grad_phi[2][1] = 1-xi;
-               grad_phi[3][1] = xi;
-
-               break;
-             }
+         case 2:
+         {
+           double xi  = q_points[i](0);
+           double eta = q_points[i](1);
+
+                                            // Define shape functions on reference element
+                                            // we consider a bi-linear mapping
+           phi[0] = (1. - xi) * (1. - eta);
+           phi[1] = xi * (1. - eta);
+           phi[2] = (1. - xi) * eta;
+           phi[3] = xi * eta;
+
+           grad_phi[0][0] = (-1. + eta);
+           grad_phi[1][0] = (1. - eta);
+           grad_phi[2][0] = -eta;
+           grad_phi[3][0] = eta;
+
+           grad_phi[0][1] = (-1. + xi);
+           grad_phi[1][1] = -xi;
+           grad_phi[2][1] = 1-xi;
+           grad_phi[3][1] = xi;
 
-             default:
-                   Assert (false, ExcNotImplemented());
-           }
+           break;
+         }
+
+         default:
+               Assert (false, ExcNotImplemented());
+       }
 
 
-         Tensor<2,dim> jacobian;
+      Tensor<2,dim> jacobian;
 
-                                          // Calculate Jacobian of transformation
-         for (unsigned int d=0; d<dim; ++d)
-           for (unsigned int e=0; e<dim; ++e)
-                       {
-             for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
-                               {
+                                      // Calculate Jacobian of transformation
+      for (unsigned int d=0; d<dim; ++d)
+       for (unsigned int e=0; e<dim; ++e)
+         {
+           for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+             {
                jacobian[d][e] += grad_phi[j][e] * v[j](d);
-                               }
-                       }
-
-         double detJ = determinant(jacobian);
-         xfem_weights.push_back (W[i] * detJ);
-
-               // Map integration points from reference element to subcell of reference element
-               Point<dim> q_prime;
-               for (unsigned int d=0; d<dim; ++d)
-                 for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
-                   q_prime[d] += v[j](d) * phi[j];
-               xfem_points.push_back(q_prime);
-       }
+             }
+         }
+
+      double detJ = determinant(jacobian);
+      xfem_weights.push_back (W[i] * detJ);
+
+                                      // Map integration points from reference element to subcell of reference element
+      Point<dim> q_prime;
+      for (unsigned int d=0; d<dim; ++d)
+       for (unsigned int j = 0; j<GeometryInfo<dim>::vertices_per_cell; j++)
+         q_prime[d] += v[j](d) * phi[j];
+      xfem_points.push_back(q_prime);
+    }
 
 }
 
@@ -848,7 +976,7 @@ void LaplaceProblem<dim>::run ()
          static const HyperBallBoundary<dim> boundary;
          triangulation.set_boundary (0, boundary);
 
-         triangulation.refine_global (1);
+         triangulation.refine_global (3);
        }
       else
        refine_grid ();

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.