]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added the lines to adjust almost 1 and almost 0 entries
authorkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 31 Jul 2006 01:27:09 +0000 (01:27 +0000)
committerkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 31 Jul 2006 01:27:09 +0000 (01:27 +0000)
in the face/subface interpolation matrices to the
method for subfaces.

git-svn-id: https://svn.dealii.org/trunk@13534 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_q.cc

index cd33fd99a298c5017073708be07f2dce66e45020..24767f8841a1d494c2e5ed5fe4ae60dc5e363214 100644 (file)
@@ -493,11 +493,17 @@ get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
          ExcDimensionMismatch (interpolation_matrix.m(),
                                source_fe.dofs_per_face));
   
-                                       // generate a point on this
-                                       // cell and evaluate the
-                                       // shape functions there
+                                   // generate a point on this
+                                   // cell and evaluate the
+                                   // shape functions there
   Quadrature<dim-1> quad_face_support (source_fe.get_unit_face_support_points ());
 
+                                  // Rule of thumb for FP accuracy,
+                                  // that can be expected for a
+                                  // given polynomial degree.
+                                  // This value is used to cut
+                                  // off values close to zero.
+  double eps = 2e-14*this->degree*(dim-1);
 
                                   // compute the interpolation
                                   // matrix by simply taking the
@@ -510,14 +516,23 @@ get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
       Point<dim> p = QProjector<dim>::project_to_subface (quad_face_support, 0, subface).point (i);
 
       for (unsigned int j=0; j<this->dofs_per_face; ++j)
-       interpolation_matrix(j,i) = this->shape_value (this->face_to_cell_index(j, 0), p);
-    }
+       { 
+         double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
 
-                                   // cut off very small values
-  for (unsigned int i=0; i<this->dofs_per_face; ++i)
-    for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
-      if (std::fabs(interpolation_matrix(i,j)) < 1e-15)
-        interpolation_matrix(i,j) = 0.;
+                                          // Correct the interpolated
+                                          // value. I.e. if it is close
+                                          // to 1 or 0, make it exactly
+                                          // 1 or 0. Unfortunately, this
+                                          // is required to avoid problems
+                                          // with higher order elements.
+         if (fabs (matrix_entry - 1.0) < eps)
+           matrix_entry = 1.0;
+         if (fabs (matrix_entry) < eps)
+           matrix_entry = 0.0;
+
+         interpolation_matrix(j,i) = matrix_entry;
+       }
+    }
 
                                   // make sure that the column sum of
                                   // each of the matrices is 1 at

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.