From: Wolfgang Bangerth <bangerth@math.tamu.edu> Date: Sat, 3 Aug 2013 21:28:49 +0000 (+0000) Subject: Fix bug #76 by removing an assertion that incorrectly triggered. The assertion was... X-Git-Tag: v8.1.0~1151 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a288ed5deba68bce6fd719699babad4c054ccbf9;p=dealii.git Fix bug #76 by removing an assertion that incorrectly triggered. The assertion was also pretty expensive. git-svn-id: https://svn.dealii.org/trunk@30219 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 401a70d3ae..b411a0066f 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -44,8 +44,15 @@ inconvenience this causes. <h3>Specific improvements</h3> <ol> + <li>Fixed: An assertion wrongly triggered in + DoFTools::make_hanging_node_constraints when used with a particular + combination of FESystem elements containing FE_Nothing. This is now fixed. + <br> + (Denis Davydov, Wolfgang Bangerth, 2013/08/01) + </li> + <li> - add has_ghost_elements() for PETScWrappers::MPI::BlockVector and + New: Add has_ghost_elements() for PETScWrappers::MPI::BlockVector and TrilinosWrappers::MPI::BlockVector. <br> (Timo Heister, 2013/08/01) diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index 3f7e6f3f39..bcc29284f5 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -2692,39 +2692,15 @@ namespace DoFTools // another possibility is what happens in crash_13. // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs. // FESystem(FE_Q(1),FE_DGQ(1)). neither of them - // dominates the other. the point is that it - // doesn't matter, since hopefully for this case, - // both sides' dofs should have been unified. - // - // make sure this is actually true. this actually - // only matters, of course, if either of the two - // finite elements actually do have dofs on the - // face - if ((cell->get_fe().dofs_per_face != 0) - || - (cell->neighbor(face)->get_fe().dofs_per_face != 0)) - { - Assert (cell->get_fe().dofs_per_face - == - cell->neighbor(face)->get_fe().dofs_per_face, - ExcNotImplemented()); - - // (ab)use the master and slave dofs arrays for - // a moment here - master_dofs.resize (cell->get_fe().dofs_per_face); - cell->face(face) - ->get_dof_indices (master_dofs, - cell->active_fe_index ()); - - slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face); - cell->face(face) - ->get_dof_indices (slave_dofs, - cell->neighbor(face)->active_fe_index ()); - - for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i) - AssertDimension (master_dofs[i], slave_dofs[i]); - } - + // dominates the other. + // + // a final possibility is that we have something like + // FESystem(FE_Q(1),FE_Q(1)) vs + // FESystem(FE_Q(1),FE_Nothing()), see + // hp/fe_nothing_18/19. + // + // in any case, the point is that it doesn't + // matter. there is nothing to do here. break; } diff --git a/tests/hp/fe_nothing_18.cc b/tests/hp/fe_nothing_18.cc new file mode 100644 index 0000000000..3bd0683d5d --- /dev/null +++ b/tests/hp/fe_nothing_18.cc @@ -0,0 +1,383 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test by Denis Davydov, see bug report 76: can't create hanging node +// constraints for a particular combination of finite elements + +#include "../tests.h" +#include <deal.II/base/function.h> +#include <deal.II/base/logstream.h> + +#include <deal.II/base/quadrature_lib.h> //needed for assembling the matrix using quadrature on each cell +#include <deal.II/base/function.h> +#include <deal.II/base/logstream.h> +#include <deal.II/lac/block_sparse_matrix.h> +#include <deal.II/lac/block_vector.h> +#include <deal.II/lac/vector.h> +#include <deal.II/lac/full_matrix.h> +#include <deal.II/lac/sparse_direct.h> //direct solver +#include <deal.II/lac/sparse_matrix.h> //definition of the sparse matrix +#include <deal.II/lac/compressed_sparsity_pattern.h> //for the intermediate sparsity pattern structure +#include <deal.II/lac/solver_cg.h> //CG solver +#include <deal.II/lac/precondition.h> //and a preconditioner +#include <deal.II/lac/constraint_matrix.h> //conform hanging nodes DoF to certain constrains to make the solution continuous +#include <deal.II/lac/iterative_inverse.h> //need for Schur matrix +#include <deal.II/grid/tria.h> //triangulation class +#include <deal.II/grid/grid_generator.h> //standard functions to generate grid +#include <deal.II/grid/grid_refinement.h> +#include <deal.II/grid/grid_tools.h> +#include <deal.II/grid/tria_accessor.h> //these two to loop over cells/faces +#include <deal.II/grid/tria_iterator.h> // ^ +#include <deal.II/grid/tria_boundary_lib.h>//to use not straightline boundaries +#include <deal.II/dofs/dof_handler.h> //associate DoF to cells/vertices/lines +#include <deal.II/dofs/dof_accessor.h> //provides information about the degrees of freedom local to a cell +#include <deal.II/dofs/dof_tools.h> //needed for the creation of sparsity patterns of sparse matrices +#include <deal.II/dofs/dof_renumbering.h> //use to renumber DoF to have a better sparsity pattern +#include <deal.II/fe/fe_values.h> //needed for assembling the matrix using quadrature on each cell +#include <deal.II/fe/fe_nothing.h> +#include <deal.II/numerics/vector_tools.h> +#include <deal.II/numerics/matrix_tools.h> +#include <deal.II/numerics/data_out.h> +#include <deal.II/numerics/error_estimator.h> +#include <deal.II/numerics/fe_field_function.h> //for evaluating solution at a given point +#include <deal.II/fe/fe_system.h> /*for vector-valued FE; vector-valude FE will be composed from the regular Q1 included below*/ +#include <deal.II/fe/fe_q.h> /*Q1 FE: 1 DoF in each vertex, none on sides...)*/ +#include <deal.II/hp/mapping_collection.h> +#include <deal.II/hp/dof_handler.h> +#include <deal.II/hp/fe_values.h> +#include <deal.II/hp/q_collection.h> +#include <fstream> +#include <iostream> +#include <cmath> //sqrt and fabs functions +#include <memory> //smart pointers +#include <algorithm> //min,max + + +template <int dim> +class ElasticProblem +{ + public: + ElasticProblem(); + ~ElasticProblem (); + void run(); + private: + enum { + u_block = 0, + lambda_block = 1 + }; + + void make_grid(); + void setup_system (); + + const unsigned int n_blocks;//total number of blocks + const unsigned int n_components;// total number of components + const unsigned int first_u_comp;//where displacements starts + const unsigned int first_lamda_comp;//where lagrage multipliers start + const unsigned int degree;//of shape functions + + std::vector<types::global_dof_index> dofs_per_block; + + Triangulation<dim> triangulation; //a triangulation object of the "dim"-dimensional domain; + hp::DoFHandler<dim> dof_handler; //is assotiated with triangulation; + + FESystem<dim> elasticity_fe; + FESystem<dim> elasticity_w_lagrange_fe; + hp::FECollection<dim> fe; //used to stack several other elements together to form one vector-valued finite element. + + ConstraintMatrix hanging_node_constraints; //object to hold hanging node constraints after refinement + ConstraintMatrix constraints; + + BlockSparsityPattern sparsity_pattern; //store sparsity pattern + + SymmetricTensor<4,dim> stress_strain_tensor; // elastic constitutive equations + + hp::QCollection<dim> quadrature_formula;//a Gauss quadrature formula to be used for integral evaluation + + const FEValuesExtractors::Vector u_fe; + const FEValuesExtractors::Vector lamda_fe; + + unsigned int id_of_lagrange_mult; + double beta1; + + SparsityPattern sparsity_pattern_nb;//needed when switching between block system and non-block. +}; + +/// +template <int dim> +SymmetricTensor<4,dim> +get_stress_strain_tensor (const double lambda, const double mu) +{ + SymmetricTensor<4,dim> tmp; + for (unsigned int i=0; i<dim; ++i) + for (unsigned int j=0; j<dim; ++j) + for (unsigned int k=0; k<dim; ++k) + for (unsigned int l=0; l<dim; ++l) + tmp[i][j][k][l] = ( + ((i==k) && (j==l) ? mu : 0.0) + + ((i==l) && (j==k) ? mu : 0.0) + + ((i==j) && (k==l) ? lambda : 0.0) + ); + return tmp; +} + +/// +/*vector-valued BC function */ +template <int dim> +class BoundaryValues : public Function<dim> +{ +public: + BoundaryValues () : Function<dim>(dim) {} + + virtual void vector_value (const Point<dim> &p, + Vector<double> &values) const; + virtual double value (const Point<dim> &p, + unsigned int compoment) const; +}; + +template <int dim> +inline +void BoundaryValues<dim>::vector_value (const Point<dim> &p, + Vector<double> &values) const +{ + values(0) = -0.001; +} + +template <int dim> +inline +double BoundaryValues<dim>::value (const Point<dim> &p, + unsigned int compoment) const +{ + return -0.001; +} + +/// +/* vector valued constrained displacements*/ +/*vector-valued RHS function */ +template <int dim> +class ConstrainValues : public Function<dim> +{ +public: + ConstrainValues (); //need to pass down to the base class of how many components the function consists; default - 1 + virtual void vector_value (const Point<dim> &p, + Vector<double> &values) const;//returns calculated values in the second argument; + virtual void vector_value_list (const std::vector<Point<dim> > &points, + std::vector<Vector<double> > &value_list) const;//values at several points at once + //prevent from calling virtual function "vector_value" to frequently +}; + +template <int dim> +ConstrainValues<dim>::ConstrainValues () +: +Function<dim> (dim) /*pass down to the main class the number of components of the function*/ +{} + +template <int dim> +inline +void ConstrainValues<dim>::vector_value (const Point<dim> &p, + Vector<double> &values) const +{ + Assert (values.size() == dim, + ExcDimensionMismatch (values.size(), dim));//check is the size of "values" is correct + //Assert (dim >= 2, ExcNotImplemented());//not implemented for 1d + values[0] = 0.0; +} + +template <int dim> +void ConstrainValues<dim>::vector_value_list (const std::vector<Point<dim> > &points, + std::vector<Vector<double> > &value_list) const +{ + Assert (value_list.size() == points.size(), + ExcDimensionMismatch (value_list.size(), points.size()));//check if input-output is consistent + + const unsigned int n_points = points.size();//number of points at which the function is to be evaluated + + for (unsigned int p=0; p<n_points; ++p) + ConstrainValues<dim>::vector_value (points[p], + value_list[p]); +} + + +/* constructor*/ +template <int dim> +ElasticProblem<dim>::ElasticProblem () +: +n_blocks(2), +n_components(dim*n_blocks), +first_u_comp(0), +first_lamda_comp(dim), +degree(1), +dofs_per_block(n_blocks), +dof_handler (triangulation), /*assotiate dof_handler to the triangulation */ +elasticity_fe (FE_Q<dim>(degree), dim, //use dim FE_Q of a given degree to represent displacements + FE_Nothing<dim>(), dim // zero extension of lagrange multipliers elsewhere in the domain +), +elasticity_w_lagrange_fe( + FE_Q<dim>(degree), dim, + FE_Q<dim>(degree), dim // same for lagrange multipliers +), +u_fe(first_u_comp), +lamda_fe(first_lamda_comp), +id_of_lagrange_mult(1), +beta1(1.0) +{ + fe.push_back (elasticity_fe);//FE index 0 + fe.push_back (elasticity_w_lagrange_fe);//FE index 1 + + quadrature_formula.push_back(QGauss<dim>(degree+2)); + quadrature_formula.push_back(QGauss<dim>(degree+2)); + + stress_strain_tensor + = get_stress_strain_tensor<dim> (/*lambda = */ 9.695e10, + /*mu = */ 7.617e10); +} + +template <int dim> +ElasticProblem<dim>::~ElasticProblem () +{ + dof_handler.clear (); +} + + +template <int dim> +void ElasticProblem<dim>::make_grid () +{ + + Triangulation<dim> triangulationL; + Triangulation<dim> triangulationR; + GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,1]^d domain + GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,1]^d domain + Point<dim> shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulationR); + + const unsigned int n_faces_per_cell = GeometryInfo< dim >::faces_per_cell; + + triangulationL.begin_active()->set_material_id(0); + triangulationR.begin_active()->set_material_id(id_of_lagrange_mult); + + for (unsigned int i=0;i<n_faces_per_cell;i++) { + triangulationL.begin_active()->face(i)->set_boundary_indicator(i); + triangulationR.begin_active()->face(i)->set_boundary_indicator(n_faces_per_cell+i); + } + + GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation); + + //triangulation.refine_global (2); //refine twice globally before solving +} + + +template <int dim> +void ElasticProblem<dim>::setup_system () +{ + std::vector<unsigned int> block_component(n_components, u_block);// init to represent u everywhere + for (unsigned int i = 0; i < dim; i++) + block_component[i+dim] = lambda_block; + + //(1) set active FE indices based in material id... + typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + unsigned int n_lagrange_cells = 0; + unsigned int n_elasticity_cells = 0; + for (; cell!=endc; ++cell) //loop over all cells + { + if (cell->material_id() == id_of_lagrange_mult ) { + cell->set_active_fe_index(1); // for lagrangian cells.. + n_lagrange_cells++; + } else { + cell->set_active_fe_index(0); + n_elasticity_cells++; + } + } + deallog<<" number of cells (L/E): "<< n_lagrange_cells << "; "<<n_elasticity_cells <<std::endl; + Assert (n_lagrange_cells > 0, ExcInternalError()); //there should be at least 1 cell! Otherwise DoFHanlder crashes with 0 dofs for block 2! + // + //(2) distribute DoFs + dof_handler.distribute_dofs (fe); /*enumerate DoF based on fe object (which knows about shape-functions used, their order and dimension*/ + DoFRenumbering::Cuthill_McKee (dof_handler);//do renumberring; must be done right after distributing DoF !!! + DoFRenumbering::component_wise(dof_handler,block_component); + DoFTools::count_dofs_per_block(dof_handler,dofs_per_block,block_component); + + deallog<<"dofs per block: U="<<dofs_per_block[u_block]<<" L="<<dofs_per_block[lambda_block]<<std::endl; + + //related to hanging nodes and refinement: + //(3) create hanging nodes constrains (also add normal diricle BC here? + hanging_node_constraints.clear ();//clear all content that may left from previous calculations + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints);//fill the object representing constraints + hanging_node_constraints.close ();//close the object - sort and rearrange constraints for effective calculations + + //(4) create block sparsity pattern (define which elements in sparse matrix are non-zero; prescribe coupling between blocks) + // following step-22 use of simple compressed block sparsity pattern for efficiency + BlockCompressedSparsityPattern compressed_sparsity_pattern(n_blocks, + n_blocks); + + compressed_sparsity_pattern.block(u_block ,u_block ).reinit(dofs_per_block[u_block] ,dofs_per_block[u_block]); + compressed_sparsity_pattern.block(u_block ,lambda_block).reinit(dofs_per_block[u_block] ,dofs_per_block[lambda_block]); + compressed_sparsity_pattern.block(lambda_block,u_block ).reinit(dofs_per_block[lambda_block],dofs_per_block[u_block]); + compressed_sparsity_pattern.block(lambda_block,lambda_block).reinit(dofs_per_block[lambda_block],dofs_per_block[lambda_block]); + + compressed_sparsity_pattern.collect_sizes(); + + Table<2, DoFTools::Coupling> coupling(n_components, n_components); + for (unsigned int ii = 0; ii< n_components;++ii) + for (unsigned int jj=0; jj< n_components; ++jj) + { + if ( (block_component[ii]==lambda_block) && + (block_component[jj]==lambda_block) + ) + coupling[ii][jj] = DoFTools::none;//diagonal = 0 + else + coupling[ii][jj] = DoFTools::always;//full coupling (u,u), (u,lambda) + } + + hanging_node_constraints.condense (compressed_sparsity_pattern); + + DoFTools::make_sparsity_pattern (dof_handler, + coupling, + compressed_sparsity_pattern, + hanging_node_constraints, + false); + + compressed_sparsity_pattern.print (deallog.get_file_stream()); + hanging_node_constraints.print (deallog.get_file_stream()); +} + + + +template <int dim> +void ElasticProblem<dim>::run () +{ + make_grid (); + setup_system (); +} + + +int main() +{ + std::ofstream logfile("fe_nothing_18/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + ElasticProblem <2> elastic_problem; + elastic_problem.run(); + + deallog << "OK" << std::endl; +} + + + diff --git a/tests/hp/fe_nothing_18/cmp/generic b/tests/hp/fe_nothing_18/cmp/generic new file mode 100644 index 0000000000..d443ae5369 --- /dev/null +++ b/tests/hp/fe_nothing_18/cmp/generic @@ -0,0 +1,24 @@ + +DEAL:: number of cells (L/E): 1; 1 +DEAL::dofs per block: U=12 L=8 +[0,0,1,2,3,4,5,6,7] +[1,0,1,2,3,4,5,6,7] +[2,0,1,2,3,4,5,6,7] +[3,0,1,2,3,4,5,6,7] +[4,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[5,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[6,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[7,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[8,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[9,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[10,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[11,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] +[12,4,5,6,7,8,9,10,11] +[13,4,5,6,7,8,9,10,11] +[14,4,5,6,7,8,9,10,11] +[15,4,5,6,7,8,9,10,11] +[16,4,5,6,7,8,9,10,11] +[17,4,5,6,7,8,9,10,11] +[18,4,5,6,7,8,9,10,11] +[19,4,5,6,7,8,9,10,11] +DEAL::OK diff --git a/tests/hp/fe_nothing_19.cc b/tests/hp/fe_nothing_19.cc new file mode 100644 index 0000000000..5ac74f5a69 --- /dev/null +++ b/tests/hp/fe_nothing_19.cc @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// a redux of fe_nothing_18 + +#include "../tests.h" +#include <deal.II/base/function.h> +#include <deal.II/base/logstream.h> + +#include <deal.II/grid/tria.h> //triangulation class +#include <deal.II/grid/grid_generator.h> //standard functions to generate grid +#include <deal.II/dofs/dof_accessor.h> +#include <deal.II/dofs/dof_tools.h> +#include <deal.II/fe/fe_nothing.h> +#include <deal.II/fe/fe_system.h> +#include <deal.II/fe/fe_q.h> +#include <deal.II/hp/dof_handler.h> +#include <deal.II/hp/q_collection.h> +#include <fstream> +#include <iostream> + +template <int dim> +void test() +{ + Triangulation<dim> triangulation; + GridGenerator::hyper_cube (triangulation, -1,0); + triangulation.refine_global (1); + + FESystem<dim> elasticity_fe (FE_Q<dim>(1), 1, + FE_Nothing<dim>(), 1); + FESystem<dim> elasticity_w_lagrange_fe(FE_Q<dim>(1), 1, + FE_Q<dim>(1), 1); + hp::FECollection<dim> fe; + fe.push_back (elasticity_fe); + fe.push_back (elasticity_w_lagrange_fe); + + hp::DoFHandler<dim> dof_handler (triangulation); + typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(); + cell->set_active_fe_index(0); + (++cell)->set_active_fe_index(1); + + dof_handler.distribute_dofs (fe); + ConstraintMatrix hanging_node_constraints; + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + + // print constraints. there shouldn't be any + deallog << hanging_node_constraints.n_constraints() << std::endl; + hanging_node_constraints.print (deallog.get_file_stream()); +} + + +int main() +{ + std::ofstream logfile("fe_nothing_19/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1>(); + test<2>(); + test<3>(); +} + + + diff --git a/tests/hp/fe_nothing_19/cmp/generic b/tests/hp/fe_nothing_19/cmp/generic new file mode 100644 index 0000000000..70626d7e3f --- /dev/null +++ b/tests/hp/fe_nothing_19/cmp/generic @@ -0,0 +1,4 @@ + +DEAL::0 +DEAL::0 +DEAL::0