]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add projection to boundary.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 22 May 1998 07:40:04 +0000 (07:40 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 22 May 1998 07:40:04 +0000 (07:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@331 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc

index 69c917d04e6c945cdff8b0890c56d02d7552480f..6880086717f402a3c1b1817df736916f60021444 100644 (file)
@@ -107,14 +107,8 @@ void Assembler<dim>::assemble (const Equation<dim> &equation) {
                    boundary);
   const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs;
 
-                                  // clear cell matrix
   if (assemble_matrix)
-    for (unsigned int i=0; i<n_dofs; ++i)
-      for (unsigned int j=0; j<n_dofs; ++j)
-       cell_matrix(i,j) = 0;
-  
-
-                                  // clear cell vector
+    cell_matrix.clear ();
   if (assemble_rhs)
     cell_vector.clear ();
   
index 36fc572414677cccad8aca3a7ceda390ecec9853..2be140e98f4eb2785a568dbde76c2035a31af5cd 100644 (file)
@@ -4,6 +4,7 @@
 #include <grid/dof.h>
 #include <grid/dof_accessor.h>
 #include <grid/tria_iterator.h>
+#include <grid/geometry_info.h>
 #include <fe/quadrature.h>
 #include <fe/fe.h>
 #include <fe/fe_values.h>
@@ -11,6 +12,9 @@
 #include <numerics/assembler.h>
 #include <lac/dsmatrix.h>
 
+#include <algorithm>
+#include <set>
+
 
 
 template <int dim>
@@ -74,6 +78,20 @@ void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim>    &dof,
 
 
 
+void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1>    &,
+                                                   const FiniteElement<1> &,
+                                                   const Quadrature<0>    &,
+                                                   const Boundary<1>      &,
+                                                   dSMatrix               &,
+                                                   const FunctionMap      &,
+                                                   dVector                &,
+                                                   vector<int>            &,
+                                                   const Function<1>      *) {
+  Assert (false, ExcNotImplemented());
+};
+
+
+
 template <int dim>
 void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &dof,
                                                      const FiniteElement<dim> &fe,
@@ -82,8 +100,141 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
                                                      dSMatrix                 &matrix,
                                                      const FunctionMap        &rhs,
                                                      dVector                  &rhs_vector,
-                                                     vector<int>              &vec_to_dof_mapping,
-                                                     const Function<dim>      *a) {};
+                                                     vector<int>              &dof_to_boundary_mapping,
+                                                     const Function<dim>      *a) {
+  Assert (matrix.n() == dof.n_boundary_dofs(rhs), ExcInternalError());
+  Assert (matrix.n() == matrix.m(), ExcInternalError());
+  Assert (matrix.n() == rhs_vector.size(), ExcInternalError());
+  Assert (rhs.size() != 0, ExcInternalError());
+  Assert (dof.get_selected_fe() == fe, ExcInternalError());
+  Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), ExcInternalError());
+  Assert (*max_element(dof_to_boundary_mapping.begin(),dof_to_boundary_mapping.end()) ==
+         (signed int)matrix.n()-1,
+         ExcInternalError());
+  
+  const unsigned int dofs_per_cell = fe.total_dofs,
+                    dofs_per_face = fe.dofs_per_face;
+  dFMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
+  dVector  cell_vector(dofs_per_cell);
+  
+  
+  UpdateFlags update_flags = UpdateFlags (update_JxW_values | update_q_points);
+  FEFaceValues<dim> fe_values (fe, q, update_flags);
+  
+  DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
+                                       endc = dof.end ();
+  for (; cell!=endc; ++cell)
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+                                      // check if this face is on that part of
+                                      // the boundary we are interested in
+      if (rhs.find(cell->face(face)->boundary_indicator()) != rhs.end())
+       {
+         cell_matrix.clear ();
+         cell_vector.clear ();
+         
+         fe_values.reinit (cell, face, fe, boundary);
+
+         const dFMatrix       &values    = fe_values.get_shape_values ();
+         const vector<double> &weights   = fe_values.get_JxW_values ();
+         vector<double>        rhs_values (fe_values.n_quadrature_points);
+         rhs.find(cell->face(face)->boundary_indicator())
+           ->second->value_list (fe_values.get_quadrature_points(), rhs_values);
+         
+         if (a != 0)
+           {
+             vector<double> coefficient_values (fe_values.n_quadrature_points);
+             a->value_list (fe_values.get_quadrature_points(), coefficient_values);
+             for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+               for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+                 {
+                   for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+                     cell_matrix(i,j) += (values(i,point) *
+                                          values(j,point) *
+                                          weights[point] *
+                                          coefficient_values[point]);
+                   cell_vector(i) += values(i,point) *
+                                     rhs_values[point] *
+                                     weights[point];
+                 };
+           }
+         else
+           for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+             for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+               {
+                 for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+                   cell_matrix(i,j) += (values(i,point) *
+                                        values(j,point) *
+                                        weights[point]);
+                 cell_vector(i) += values(i,point) *
+                                   rhs_values[point] *
+                                   weights[point];
+               };
+
+
+                                          // now transfer cell matrix and vector
+                                          // to the whole boundary matrix
+                                          //
+                                          // in the following: dof[i] holds the
+                                          // global index of the i-th degree of
+                                          // freedom on the present cell. If it
+                                          // is also a dof on the boundary, it
+                                          // must be a nonzero entry in the
+                                          // dof_to_boundary_mapping and then
+                                          // the boundary index of this dof is
+                                          // dof_to_boundary_mapping[dof[i]].
+                                          //
+                                          // if dof[i] is not on the boundary,
+                                          // it should be zero on the boundary
+                                          // therefore on all quadrature
+                                          // points and finally all of its
+                                          // entries in the cell matrix and
+                                          // vector should be zero. If not, we
+                                          // throw an error
+                                          //
+                                          // the main problem here is that the
+                                          // matrix or vector entry should also
+                                          // be zero if the degree of freedom
+                                          // dof[i] is on the boundary, but not
+                                          // on the present face, i.e. on
+                                          // another face of the same cell also
+                                          // on the boundary. We can therefore
+                                          // not rely on the
+                                          // dof_to_boundary_mapping[dof[i]]
+                                          // being !=-1, we really have to
+                                          // determine whether dof[i] is a
+                                          // dof on the present face. We do so
+                                          // by getting the dofs on the
+                                          // face into #dofs_on_face_vector#,
+                                          // a vector as always. Usually,
+                                          // searching in a vector is
+                                          // inefficient, so we copy the dofs
+                                          // into a set, which enables binary
+                                          // searches.
+         vector<int> dofs (dofs_per_cell);
+         cell->get_dof_indices (dofs);
+
+         vector<int> dofs_on_face_vector (dofs_per_face);
+         cell->face(face)->get_dof_indices (dofs_on_face_vector);
+         set<int> dofs_on_face (dofs_on_face_vector.begin(),
+                                dofs_on_face_vector.end());
+         
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if ((dofs_on_face.find(dofs[i]) != dofs_on_face.end()) &&
+                 (dofs_on_face.find(dofs[j]) != dofs_on_face.end()))
+               matrix.add(dof_to_boundary_mapping[dofs[i]],
+                          dof_to_boundary_mapping[dofs[j]],
+                          cell_matrix(i,j));
+             else
+               Assert (cell_matrix(i,j)==0, ExcInternalError ());
+         
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           if (dofs_on_face.find(dofs[j]) != dofs_on_face.end())
+             rhs_vector(dof_to_boundary_mapping[dofs[j]]) += cell_vector(j);
+           else
+             Assert (cell_vector(j)==0, ExcInternalError());
+       };
+};
 
 
 
index b395244d4313a9577ef17fe60b481060b381d035..c207759144596df12e91a0c72ceb4e0f2a808ddf 100644 (file)
@@ -242,7 +242,15 @@ VectorTools<dim>::project_boundary_values (const DoFHandler<dim>    &dof,
                                   // in this case the boundary mass matrix has
                                   // to be condensed and the solution is to
                                   // be distributed afterwards, which is not
-                                  // yet implemented
+                                  // yet implemented. The reason for this is
+                                  // that we cannot simply use the #condense#
+                                  // family of functions, since the matrices
+                                  // and vectors do not use the global
+                                  // numbering but rather the boundary
+                                  // numbering, i.e. the condense function
+                                  // needs to use another indirection. There
+                                  // should be not many technical problems,
+                                  // but it needs to be implemented
   if (dim<3)
     sparsity.compress();
   else
@@ -257,7 +265,34 @@ VectorTools<dim>::project_boundary_values (const DoFHandler<dim>    &dof,
   MatrixTools<dim>::create_boundary_mass_matrix (dof, fe, q, boundary,
                                                 mass_matrix, boundary_functions,
                                                 rhs, dof_to_boundary_mapping);
+
+                                  // same thing as above: if dim>=3 we need
+                                  // to consider constraints
+  Assert (dim<3, ExcNotImplemented());
+
   
+  dVector boundary_projection (rhs.size());
+
+  int    max_iter  = 1000;
+  double tolerance = 1.e-16;
+  Control                          control1(max_iter,tolerance);
+  PrimitiveVectorMemory<dVector>   memory(rhs.size());
+  CG<dSMatrix,dVector>             cg(control1,memory);
+
+                                  // solve
+  cg (mass_matrix, boundary_projection, rhs);
+
+                                  // fill in boundary values
+  for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
+    if (dof_to_boundary_mapping[i] != -1)
+                                      // this dof is on one of the
+                                      // interesting boundary parts
+                                      //
+                                      // remember: #i# is the global dof
+                                      // number, #dof_to_boundary_mapping[i]#
+                                      // is the number on the boundary and
+                                      // thus in the solution vector
+      boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
 };
 
 

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.