]> https://gitweb.dealii.org/ - dealii.git/commitdiff
allow interpolate for combination FENothing/FE_DGQ 16496/head
authorMagdalena Schreter <magdalena.schreter@tum.de>
Thu, 18 Jan 2024 21:15:10 +0000 (22:15 +0100)
committerMagdalena Schreter <magdalena.schreter@tum.de>
Thu, 18 Jan 2024 21:49:57 +0000 (22:49 +0100)
doc/news/changes/minor/20240118Schreter [new file with mode: 0644]
source/fe/fe_dgq.cc
tests/fe/dgq_1.cc
tests/fe/dgq_1.output
tests/hp/solution_transfer_03.cc

diff --git a/doc/news/changes/minor/20240118Schreter b/doc/news/changes/minor/20240118Schreter
new file mode 100644 (file)
index 0000000..3021b47
--- /dev/null
@@ -0,0 +1,4 @@
+Changed:
+FE_DGQ::get_interpolation_matrix() now also takes as a source element FENothing.
+<br>
+(Magdalena Schreter, Peter Munch, 2024/01/18)
index 3f24eb387584e0e2ebe12ddfd9222c9fb8418351..f1d15b8667db4a5843489daae1f50b131d1e0249 100644 (file)
@@ -274,78 +274,93 @@ FE_DGQ<dim, spacedim>::get_interpolation_matrix(
   const FiniteElement<dim, spacedim> &x_source_fe,
   FullMatrix<double>                 &interpolation_matrix) const
 {
-  // this is only implemented, if the
-  // source FE is also a
-  // DGQ element
-  using FE = FiniteElement<dim, spacedim>;
-  AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
-               nullptr),
-              typename FE::ExcInterpolationNotImplemented());
-
-  // ok, source is a Q element, so
-  // we will be able to do the work
-  const FE_DGQ<dim, spacedim> &source_fe =
-    dynamic_cast<const FE_DGQ<dim, spacedim> &>(x_source_fe);
-
-  Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
-         ExcDimensionMismatch(interpolation_matrix.m(),
-                              this->n_dofs_per_cell()));
-  Assert(interpolation_matrix.n() == source_fe.n_dofs_per_cell(),
-         ExcDimensionMismatch(interpolation_matrix.n(),
-                              source_fe.n_dofs_per_cell()));
-
-
-  // compute the interpolation
-  // matrices in much the same way as
-  // we do for the embedding matrices
-  // from mother to child.
-  FullMatrix<double> cell_interpolation(this->n_dofs_per_cell(),
-                                        this->n_dofs_per_cell());
-  FullMatrix<double> source_interpolation(this->n_dofs_per_cell(),
-                                          source_fe.n_dofs_per_cell());
-  FullMatrix<double> tmp(this->n_dofs_per_cell(), source_fe.n_dofs_per_cell());
-  for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
+  // go through the list of elements we can interpolate from
+  if (const FE_DGQ<dim, spacedim> *source_fe =
+        dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe))
     {
-      // generate a point on this
-      // cell and evaluate the
-      // shape functions there
-      const Point<dim> p = this->unit_support_points[j];
-      for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
-        cell_interpolation(j, i) = this->poly_space->compute_value(i, p);
-
-      for (unsigned int i = 0; i < source_fe.n_dofs_per_cell(); ++i)
-        source_interpolation(j, i) = source_fe.poly_space->compute_value(i, p);
-    }
+      Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
+             ExcDimensionMismatch(interpolation_matrix.m(),
+                                  this->n_dofs_per_cell()));
+      Assert(interpolation_matrix.n() == source_fe->n_dofs_per_cell(),
+             ExcDimensionMismatch(interpolation_matrix.n(),
+                                  source_fe->n_dofs_per_cell()));
+
+
+      // compute the interpolation
+      // matrices in much the same way as
+      // we do for the embedding matrices
+      // from mother to child.
+      FullMatrix<double> cell_interpolation(this->n_dofs_per_cell(),
+                                            this->n_dofs_per_cell());
+      FullMatrix<double> source_interpolation(this->n_dofs_per_cell(),
+                                              source_fe->n_dofs_per_cell());
+      FullMatrix<double> tmp(this->n_dofs_per_cell(),
+                             source_fe->n_dofs_per_cell());
+      for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
+        {
+          // generate a point on this
+          // cell and evaluate the
+          // shape functions there
+          const Point<dim> p = this->unit_support_points[j];
+          for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
+            cell_interpolation(j, i) = this->poly_space->compute_value(i, p);
+
+          for (unsigned int i = 0; i < source_fe->n_dofs_per_cell(); ++i)
+            source_interpolation(j, i) =
+              source_fe->poly_space->compute_value(i, p);
+        }
 
-  // then compute the
-  // interpolation matrix matrix
-  // for this coordinate
-  // direction
-  cell_interpolation.gauss_jordan();
-  cell_interpolation.mmult(interpolation_matrix, source_interpolation);
+      // then compute the
+      // interpolation matrix matrix
+      // for this coordinate
+      // direction
+      cell_interpolation.gauss_jordan();
+      cell_interpolation.mmult(interpolation_matrix, source_interpolation);
 
-  // cut off very small values
-  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
-    for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j)
-      if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
-        interpolation_matrix(i, j) = 0.;
+      // cut off very small values
+      for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
+        for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
+          if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
+            interpolation_matrix(i, j) = 0.;
 
 #ifdef DEBUG
-  // make sure that the row sum of
-  // each of the matrices is 1 at
-  // this point. this must be so
-  // since the shape functions sum up
-  // to 1
-  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
-    {
-      double sum = 0.;
-      for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j)
-        sum += interpolation_matrix(i, j);
+      // make sure that the row sum of
+      // each of the matrices is 1 at
+      // this point. this must be so
+      // since the shape functions sum up
+      // to 1
+      for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
+        {
+          double sum = 0.;
+          for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
+            sum += interpolation_matrix(i, j);
 
-      Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim,
-             ExcInternalError());
-    }
+          Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim,
+                 ExcInternalError());
+        }
 #endif
+    }
+  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe))
+    {
+      // the element we want to interpolate from is an FE_Nothing. this
+      // element represents a function that is constant zero and has no
+      // degrees of freedom, so the interpolation is simply a multiplication
+      // with a n_dofs x 0 matrix. there is nothing to do here
+
+      // we would like to verify that the number of rows and columns of
+      // the matrix equals this->n_dofs_per_cell() and zero. unfortunately,
+      // whenever we do FullMatrix::reinit(m,0), it sets both rows and
+      // columns to zero, instead of m and zero. thus, only test the
+      // number of columns
+      Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
+             ExcDimensionMismatch(interpolation_matrix.m(),
+                                  x_source_fe.n_dofs_per_cell()));
+    }
+  else
+    AssertThrow(
+      false,
+      (typename FiniteElement<dim,
+                              spacedim>::ExcInterpolationNotImplemented()));
 }
 
 
index 879f0967c1504b1493bc28cd8f67133f08260671..974ad0fea171c76432bc7e1352fd8366c84f28a5 100644 (file)
@@ -21,6 +21,7 @@
 // result doesn't change
 
 #include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_tools.h>
 
 #include <string>
@@ -56,6 +57,31 @@ test(const unsigned int degree1, const unsigned int degree2)
 }
 
 
+template <int dim>
+void
+test_fe_nothing()
+{
+  deallog << "FE_DGQ<" << dim << ">(1)"
+          << " to FE_Nothing<" << dim << ">()" << std::endl;
+
+  FE_DGQ<dim>     fe1(1);
+  FE_Nothing<dim> fe2;
+
+  FullMatrix<double> m;
+  fe1.get_interpolation_matrix(fe2, m);
+
+  for (unsigned int i = 0; i < m.m(); ++i)
+    {
+      for (unsigned int j = 0; j < m.n(); ++j)
+        deallog << m(i, j) << ' ';
+
+      deallog << std::endl;
+    }
+
+  deallog << std::endl;
+}
+
+
 int
 main()
 {
@@ -74,5 +100,9 @@ main()
     for (unsigned int degree2 = 0; degree2 <= 2; ++degree2)
       test<3>(degree1, degree2);
 
+  test_fe_nothing<1>();
+  test_fe_nothing<2>();
+  test_fe_nothing<3>();
+
   return 0;
 }
index 8686f9a7111dd1b34ff8f926e574f4661d4e8550..1b6ca44b963a641fc93b43fc0c7810cfd75e42c7 100644 (file)
@@ -402,3 +402,9 @@ DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.
 DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 
 DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 
 DEAL::
+DEAL::FE_DGQ<1>(1) to FE_Nothing<1>()
+DEAL::
+DEAL::FE_DGQ<2>(1) to FE_Nothing<2>()
+DEAL::
+DEAL::FE_DGQ<3>(1) to FE_Nothing<3>()
+DEAL::
index f2f1c52ae674fe59f117f724cabc3547372c1933..d4e94555b4057ff41bd56dc6414a154b2c92eb54 100644 (file)
@@ -91,7 +91,7 @@ main()
   data_out.write_vtu(deallog.get_file_stream());
 
 
-  // Interpoalte solution
+  // Interpolate solution
   SolutionTransfer<2, Vector<double>> solultion_trans(dof_handler);
   solultion_trans.prepare_for_coarsening_and_refinement(solution);
 

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.