]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix a bug where we try to access elements in a matrix that are beyond the end.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 4 Sep 2012 13:46:12 +0000 (13:46 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 4 Sep 2012 13:46:12 +0000 (13:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@26233 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/base/thread_management.h
deal.II/include/deal.II/lac/constraint_matrix.templates.h
tests/deal.II/constraints_block_01.cc [new file with mode: 0644]
tests/deal.II/constraints_block_01/cmp/generic [new file with mode: 0644]

index 56a6582e3439436e0fb000cafcca63411200eed8..62ef82026109360c83e0adae3b8d140c4a0be8d1 100644 (file)
@@ -50,7 +50,11 @@ never working correctly and it is not used.
 <h3>Specific improvements</h3>
 
 <ol>
-<li> Nothing so far.
+<li> Fixed: When applying a ConstraintMatrix to a block matrix
+where the last few rows are empty, we ran into an unrelated assertion.
+This is now fixed.
+<br>
+(Jason Sheldon, Wolfgang Bangerth, 2012/09/04)
 </ol>
 
 
index 6b5f941cfcccaefa8a613d460a9ab29e7f89406e..d2f83f0333aebf111fcba8dee47b9920584fbd08 100644 (file)
@@ -3739,7 +3739,6 @@ namespace Threads
 
     template <typename> struct TaskDescriptor;
 
-
                                      /**
                                       * The task class for TBB that is
                                       * used by the TaskDescriptor
@@ -3910,7 +3909,7 @@ namespace Threads
 
 
         template <typename> friend struct TaskEntryPoint;
-        template <typename> friend struct Task;
+        friend class dealii::Threads::Task<RT>;
     };
 
 
index ea3f98ab400a69d4ff8d9fe6d26f0c6ebee70319..db31f44eb4cfd826c38fe3fce1ea67f057db8790 100644 (file)
@@ -1768,9 +1768,9 @@ namespace internals
     const unsigned int loc_row = global_rows.local_row(i);
 
 #ifdef DEBUG
-    const unsigned int * col_ptr = sparsity.n_nonzero_elements() == 0 ? 0 :
+    const unsigned int * col_ptr = sparsity.row_length(row) == 0 ? 0 :
                                    &sparsity_struct[row_start[row]];
-    number * val_ptr = sparsity.n_nonzero_elements() == 0 ? 0 :
+    number * val_ptr = sparsity.row_length(row) == 0 ? 0 :
                        &sparse_matrix->global_entry (row_start[row]);
 #else
     const unsigned int * col_ptr = &sparsity_struct[row_start[row]];
diff --git a/tests/deal.II/constraints_block_01.cc b/tests/deal.II/constraints_block_01.cc
new file mode 100644 (file)
index 0000000..e2ed848
--- /dev/null
@@ -0,0 +1,348 @@
+//----------------------------  constraints.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2012 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  constraints.cc  ---------------------------
+
+
+// test that we fixed a case where we tried to deal with some constraints in a
+// block matrix where the last few rows of one of the blocks were empty and we
+// ran into an unrelated assertion because we were accessing something beyond
+// the end of the array
+//
+// testcase by Jason Sheldon
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <base/quadrature_lib.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <lac/vector.h>
+#include <lac/sparse_direct.h>
+#include <lac/full_matrix.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/constraint_matrix.h>
+
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+
+#include <hp/dof_handler.h>
+#include <hp/q_collection.h>
+#include <hp/fe_collection.h>
+#include <hp/fe_values.h>
+
+#include <fe/fe.h>
+#include <fe/fe_tools.h>
+#include <fe/fe_q.h>
+#include <fe/fe_nothing.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+
+#include <grid/grid_generator.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria.h>
+
+#include <fstream>
+#include <cmath>
+#include <cstdlib>
+
+
+std::ofstream logfile("constraints_block_01/output");
+
+
+int main ()
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+                                  // setting up some constants
+  const unsigned int dim = 2;
+  const unsigned int solid_dim = 2*dim;
+  const unsigned int fluid_dim = dim+1;
+  const unsigned int mesh_dim  = dim;
+  const unsigned int total_dim = solid_dim + fluid_dim + mesh_dim;
+
+                                  // make the tria and domain
+  Triangulation<dim>   tria;
+
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(1);
+
+  Point<dim> real_cell_center;
+
+                                  // create a solid domain (0) in the
+                                  // lower left hand corner of a 2x2 grid
+                                  // the rest is fluid/mesh (1)
+
+  for (Triangulation<dim>::active_cell_iterator
+        cell = tria.begin_active();
+       cell != tria.end(); ++cell)
+    {
+      real_cell_center = cell->center();
+      if (real_cell_center(0) < 0.5
+         && real_cell_center(1) < 0.5 )
+       cell->set_material_id (0); // solid
+      else
+       cell->set_material_id (1); //fluid
+    }//cell
+
+                                  //create the FE spaces for the solid and the fluid/mesh
+                                  //each are padded with FE_Nothing to be equal length
+
+  std::string solid_fe_name = "FESystem[FE_Q(2)^2-FE_Q(2)^2-FE_Nothing()^2-FE_Nothing()-FE_Nothing()^2]";
+  std::string fluid_fe_name = "FESystem[FE_Nothing()^2-FE_Nothing()^2-FE_Q(2)^2-FE_Q(1)-FE_Q(2)^2]";
+
+  hp::FECollection<dim> fe_collection;
+  FiniteElement<dim> * solid_fe = FETools::get_fe_from_name<dim>(solid_fe_name);
+  FiniteElement<dim> * fluid_fe = FETools::get_fe_from_name<dim>(fluid_fe_name);
+
+  deallog << "Solid FE Space: " << solid_fe->get_name() << std::endl;
+  deallog << "Fluid/Mesh FE Space: " << fluid_fe->get_name() << std::endl;
+
+  fe_collection.push_back(*solid_fe);
+  fe_collection.push_back(*fluid_fe);
+
+  hp::DoFHandler<dim> dh(tria);
+
+  for (hp::DoFHandler<dim>::active_cell_iterator
+        cell = dh.begin_active();
+       cell != dh.end(); ++cell)
+    {
+      if (int(cell->material_id()) == 0)
+       cell->set_active_fe_index (0);
+      else
+       cell->set_active_fe_index (1);
+    }//cell
+  dh.distribute_dofs(fe_collection);
+
+  std::vector<unsigned int> block_component (total_dim,0);
+
+  for (unsigned int comp=0; comp < total_dim; comp++)
+    {
+      if ( comp < solid_dim )
+       block_component[comp] = 0;
+      else if ( comp < solid_dim+fluid_dim )
+       block_component[comp] = 1;
+      else
+       block_component[comp] = 2;
+    }//comp
+
+  std::vector<unsigned int> dofs_per_block(3,0);//3 blocks, count dofs:
+  DoFTools::count_dofs_per_component (dh, dofs_per_block, false, block_component);
+
+  DoFRenumbering::component_wise(dh, block_component);
+
+                                  //build the sparsitypattern
+
+  BlockSparsityPattern   block_sparsity_pattern;
+  {
+    BlockCompressedSimpleSparsityPattern  csp(3,3);
+    csp.block(0,0).reinit (dofs_per_block[0],dofs_per_block[0]);//solid-solid
+    csp.block(0,1).reinit (dofs_per_block[0],dofs_per_block[1]);//solid-fluid
+    csp.block(0,2).reinit (dofs_per_block[0],dofs_per_block[2]);//solid-mesh
+
+    csp.block(1,0).reinit (dofs_per_block[1],dofs_per_block[0]);//fluid-solid
+    csp.block(1,1).reinit (dofs_per_block[1],dofs_per_block[1]);//fluid-fluid
+    csp.block(1,2).reinit (dofs_per_block[1],dofs_per_block[2]);//fluid-mesh
+
+    csp.block(2,0).reinit (dofs_per_block[2],dofs_per_block[0]);//mesh-solid
+    csp.block(2,1).reinit (dofs_per_block[2],dofs_per_block[1]);//mesh-fluid
+    csp.block(2,2).reinit (dofs_per_block[2],dofs_per_block[2]);//mesh-mesh
+
+    csp.collect_sizes();
+
+                                    //enforce coupling across cells and interface
+
+    Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(), fe_collection.n_components());
+    Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(), fe_collection.n_components());
+
+    for (unsigned int c=0; c<fe_collection.n_components(); ++c)
+      {
+       for (unsigned int d=0; d<fe_collection.n_components(); ++d)
+         {
+           if (((c<solid_dim) && (d<solid_dim)) //couples solid dims with solid dims
+               ||
+               (((c>=solid_dim) && (d>=solid_dim)) //couples fluid dims and mesh dims
+                && !((c==solid_dim+dim) && (d==solid_dim+dim)))) //fluid pressure does not couple with itself
+             {
+               cell_coupling[c][d] = DoFTools::always;
+               cell_coupling[d][c] = DoFTools::always;
+             }//cell_coupling
+
+           if ((c<solid_dim) && (d>=solid_dim)) //couples solid dims with fluid and mesh dims on the interface
+             {
+               face_coupling[c][d] = DoFTools::always;
+               face_coupling[d][c] = DoFTools::always;
+             }//face_coupling
+         }//d
+      }//c
+    DoFTools::make_flux_sparsity_pattern (dh, csp, cell_coupling, face_coupling);
+    block_sparsity_pattern.copy_from (csp);
+  }
+
+                                  // build matrices and vectors
+
+  BlockSparseMatrix<double> system_matrix(block_sparsity_pattern);
+  BlockVector<double> system_rhs(dofs_per_block);
+
+                                  // set up constraints for solution and solution update
+                                  /**
+                                   * These constraints enforce that the solid displacement
+                                   * and mesh displacement are the same on the interface.
+                                   *
+                                   * They also constrain that the solid velocity
+                                   * and fluid velocity are the same on the interface
+                                   *
+                                   * The interface check is simplified for this 2x2 case
+                                   */
+
+  ConstraintMatrix constraints;
+
+  const unsigned int   dofs_per_fl_msh_face = fluid_fe->dofs_per_face;
+  const unsigned int   dofs_per_solid_face  = solid_fe->dofs_per_face;
+  std::vector<unsigned int> fl_msh_face_dof_indices (dofs_per_fl_msh_face);
+  std::vector<unsigned int> solid_face_dof_indices  (dofs_per_solid_face );
+
+  std::vector<std::pair<unsigned int, unsigned int> > solid_fluid_pairs;
+  std::vector<std::pair<unsigned int, unsigned int> > solid_mesh_pairs;
+
+  for (hp::DoFHandler<dim>::active_cell_iterator
+        cell = dh.begin_active(); cell!=dh.end(); ++cell) //loops over the cells
+    {
+      if (int(cell->material_id()) == 1)
+       {//Only loop over cells in the fluid region
+         for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+           {
+             if (!cell->at_boundary(f))
+               {
+                 bool face_is_on_interface = false;
+                                                  //checks to see if cell face neighbor is in the solid domain
+                 if (int(cell->neighbor(f)->material_id()) == 0 )
+                   face_is_on_interface = true;
+
+                 if (face_is_on_interface)
+                   {
+                     std::vector<unsigned int> solid_disp_dof, solid_vel_dof, fluid_vel_dof, mesh_disp_dof;
+
+                     cell->face(f)->get_dof_indices (solid_face_dof_indices, 0);
+                     cell->face(f)->get_dof_indices (fl_msh_face_dof_indices, 1);
+
+                     for (unsigned int i=0; i<dofs_per_solid_face; ++i)
+                       {
+                         unsigned int comp = solid_fe->face_system_to_component_index(i).first;
+                         if (comp <dim)
+                           solid_disp_dof.push_back(i);
+                         if (comp >= dim && comp < solid_dim)
+                           solid_vel_dof.push_back(i);
+                       }//i
+                     for (unsigned int i=0; i<dofs_per_fl_msh_face; ++i)
+                       {
+                         unsigned int comp = fluid_fe->face_system_to_component_index(i).first;
+                         if (comp >= solid_dim && comp < solid_dim+dim)
+                           fluid_vel_dof.push_back(i);
+                         if (comp>=solid_dim+fluid_dim)
+                           mesh_disp_dof.push_back(i);
+                       }//i
+
+                     for (unsigned int i=0; i<solid_vel_dof.size(); ++i)
+                       {//in this example solid_vel_dof.size()==fluid_vel_dof.size()
+                         constraints.add_line  (fl_msh_face_dof_indices[fluid_vel_dof[i]]);
+                         constraints.add_entry (fl_msh_face_dof_indices[fluid_vel_dof[i]],
+                                                solid_face_dof_indices[solid_vel_dof[i]], 1.0);
+
+                         solid_fluid_pairs.push_back(std::pair<unsigned int,unsigned int>
+                                                     (solid_face_dof_indices[solid_vel_dof[i]],
+                                                      fl_msh_face_dof_indices[fluid_vel_dof[i]]));
+                       }//i
+                     for (unsigned int i=0; i<solid_disp_dof.size(); ++i)
+                       {
+                         constraints.add_line  (fl_msh_face_dof_indices[mesh_disp_dof[i]]);
+                         constraints.add_entry (fl_msh_face_dof_indices[mesh_disp_dof[i]],
+                                                solid_face_dof_indices[solid_disp_dof[i]], 1.0);
+
+                         solid_mesh_pairs.push_back(std::pair<unsigned int,unsigned int>
+                                                    (solid_face_dof_indices[solid_disp_dof[i]],
+                                                     fl_msh_face_dof_indices[mesh_disp_dof[i]]));
+                       }
+                   }//at interface check
+               }//not at boundary check
+           }//face
+       }//is this in the fluid material region?
+    }//cell
+
+  constraints.close();
+
+                                  // prints out which dofs are coupled
+  deallog<<"---------------Coupled dofs---------------"<<std::endl;
+  for (unsigned int i=0; i<solid_fluid_pairs.size(); i++)
+    {
+      deallog<<"solid dof: "<< solid_fluid_pairs[i].first
+              <<", fluid dof: "<<solid_fluid_pairs[i].second<<std::endl;
+    }
+  for (unsigned int i=0; i<solid_fluid_pairs.size(); i++)
+    {
+      deallog<<"solid dof: "<<solid_mesh_pairs[i].first
+              <<", mesh dof: "<<solid_mesh_pairs[i].second<<std::endl;
+    }
+  deallog<<"------------------------------------------"<<std::endl;
+
+
+                                  //code crashes in the fluid assembly, mocked up below
+  {//fluid assembly
+    const unsigned int dofs_per_cell   = fluid_fe->dofs_per_cell;
+
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+    for (hp::DoFHandler<dim>::active_cell_iterator
+          cell = dh.begin_active();
+        cell!=dh.end();
+        ++cell) //loops over the cells
+      {
+       if (int(cell->material_id()) == 1)//are we in the fluid region?
+         {
+           cell->get_dof_indices (local_dof_indices);
+
+                                            // normally the governing equations would be added
+                                            // to the cell matrix and right hand side here
+
+           deallog<<"I am about to distribute cell: "<<cell<<std::endl;
+           deallog<<"Cell: "<<cell<<" has it's center at "<<cell->center()<<std::endl;
+
+           constraints.distribute_local_to_global(
+             cell_matrix,
+             cell_rhs,
+             local_dof_indices,
+             system_matrix,
+             system_rhs);
+
+         }//is this in the fluid material region?
+
+      }//cell
+
+  }//assemble fluid system
+
+   return 0;
+}
+
diff --git a/tests/deal.II/constraints_block_01/cmp/generic b/tests/deal.II/constraints_block_01/cmp/generic
new file mode 100644 (file)
index 0000000..e8cfce6
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL::Solid FE Space: FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(2)^2-FE_Nothing<2>()^2-FE_Nothing<2>()-FE_Nothing<2>()^2]
+DEAL::Fluid/Mesh FE Space: FESystem<2>[FE_Nothing<2>()^2-FE_Nothing<2>()^2-FE_Q<2>(2)^2-FE_Q<2>(1)-FE_Q<2>(2)^2]
+DEAL::---------------Coupled dofs---------------
+DEAL::solid dof: 6, fluid dof: 36
+DEAL::solid dof: 7, fluid dof: 37
+DEAL::solid dof: 14, fluid dof: 42
+DEAL::solid dof: 15, fluid dof: 43
+DEAL::solid dof: 22, fluid dof: 48
+DEAL::solid dof: 23, fluid dof: 49
+DEAL::solid dof: 10, fluid dof: 58
+DEAL::solid dof: 11, fluid dof: 59
+DEAL::solid dof: 14, fluid dof: 42
+DEAL::solid dof: 15, fluid dof: 43
+DEAL::solid dof: 30, fluid dof: 71
+DEAL::solid dof: 31, fluid dof: 72
+DEAL::solid dof: 4, mesh dof: 86
+DEAL::solid dof: 5, mesh dof: 87
+DEAL::solid dof: 12, mesh dof: 90
+DEAL::solid dof: 13, mesh dof: 91
+DEAL::solid dof: 20, mesh dof: 94
+DEAL::solid dof: 21, mesh dof: 95
+DEAL::solid dof: 8, mesh dof: 104
+DEAL::solid dof: 9, mesh dof: 105
+DEAL::solid dof: 12, mesh dof: 90
+DEAL::solid dof: 13, mesh dof: 91
+DEAL::solid dof: 28, mesh dof: 114
+DEAL::solid dof: 29, mesh dof: 115
+DEAL::------------------------------------------
+DEAL::I am about to distribute cell: 1.1
+DEAL::Cell: 1.1 has it's center at 0.75 0.25
+DEAL::I am about to distribute cell: 1.2
+DEAL::Cell: 1.2 has it's center at 0.25 0.75
+DEAL::I am about to distribute cell: 1.3
+DEAL::Cell: 1.3 has it's center at 0.75 0.75

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.