Merge pull request #24 from nicola-giuliani/new_quadrature_set_up

new singular quadrature routine

Author: luca-heltai

include/bem_problem.h

@@ -128,6 +128,8 @@ class BEMProblem : public ParameterAcceptor
   /// in our computations (assemble system and compute_normals-gradients).
   void reinit();
 
+  const Quadrature<dim-1> & get_singular_quadrature(const unsigned int index) const;
+
   /// This function compute a very specific case, a double node that has a
   /// dirichlet-dirichlet condition. In this case there is a constraint for
   /// the normal derivative since we want a conitnuos velocity thus a conitnuos

source/bem_problem.cc

@@ -316,6 +316,44 @@ void BEMProblem<dim>::reinit()
 
 }
 
+
+template<>
+const Quadrature<2> &BEMProblem<3>::get_singular_quadrature(const unsigned int index) const
+{
+  Assert(index < fe->dofs_per_cell,
+         ExcIndexRange(0, fe->dofs_per_cell, index));
+
+
+
+  static std::vector<Quadrature<2> > quadratures;
+  {
+    if (quadratures.size() == 0)
+      for (unsigned int i=0; i<fe->dofs_per_cell; ++i)
+        {
+          quadratures.push_back(QSplit<2> (QDuffy (singular_quadrature_order,1.),fe->get_unit_support_points()[i]));
+        }
+  }
+
+  return quadratures[index];
+
+}
+
+template<>
+const Quadrature<1> &BEMProblem<2>::get_singular_quadrature(const unsigned int index) const
+{
+  Assert(index < fe->dofs_per_cell,
+         ExcIndexRange(0, fe->dofs_per_cell, index));
+
+  static std::vector<Quadrature<1> > quadratures;
+  if (quadratures.size() == 0)
+    for (unsigned int i=0; i<fe->dofs_per_cell; ++i)
+      {
+        quadratures.push_back(QTelles<1>(singular_quadrature_order,
+                                         fe->get_unit_support_points()[i]));
+      }
+  return quadratures[index];
+}
+
 template <int dim>
 void BEMProblem<dim>::declare_parameters (ParameterHandler &prm)
 {
@@ -566,25 +604,6 @@ void BEMProblem<dim>::assemble_system()
   dirichlet_matrix = 0;
 
 
-  std::vector<Quadrature<dim-1> > sing_quadratures;
-  for (unsigned int i=0; i<fe->dofs_per_cell; ++i)
-    {
-      if (fe->degree > 1)
-        {
-          sing_quadratures.push_back(QIterated<dim-1>(QGauss<1> (singular_quadrature_order),fe->degree));
-        }
-      else
-        {
-          sing_quadratures.push_back
-          (QTelles<dim-1>(singular_quadrature_order,
-                          fe->get_unit_support_points()[i]));
-          //
-          // Usage of alternative singular quadrature formula
-          // (QGaussOneOverR<dim-1>(singular_quadrature_order,
-          //                 fe->get_unit_support_points()[i],true));
-        }
-    }
-
 
   // Next, we initialize an FEValues
   // object with the quadrature
@@ -945,8 +964,7 @@ void BEMProblem<dim>::assemble_system()
 
                   const Quadrature<dim-1> *
                   singular_quadrature
-                    = dynamic_cast<Quadrature<dim-1>*>(
-                        &sing_quadratures[singular_index]);
+                    = &(get_singular_quadrature(singular_index));
                   Assert(singular_quadrature, ExcInternalError());
 
                   FEValues<dim-1,dim> fe_v_singular (*mapping, *fe, *singular_quadrature,