From 12a61a90b0500691d986fe1812ef299e5ecf5ebb Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 3 Aug 2013 21:28:49 +0000 Subject: [PATCH] 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 --- deal.II/doc/news/changes.h | 9 +- deal.II/source/dofs/dof_tools.cc | 42 +--- tests/hp/fe_nothing_18.cc | 383 +++++++++++++++++++++++++++++ tests/hp/fe_nothing_18/cmp/generic | 24 ++ tests/hp/fe_nothing_19.cc | 81 ++++++ tests/hp/fe_nothing_19/cmp/generic | 4 + 6 files changed, 509 insertions(+), 34 deletions(-) create mode 100644 tests/hp/fe_nothing_18.cc create mode 100644 tests/hp/fe_nothing_18/cmp/generic create mode 100644 tests/hp/fe_nothing_19.cc create mode 100644 tests/hp/fe_nothing_19/cmp/generic 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.

Specific improvements

    +
  1. 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. +
    + (Denis Davydov, Wolfgang Bangerth, 2013/08/01) +
    2. +
  2. - add has_ghost_elements() for PETScWrappers::MPI::BlockVector and + New: Add has_ghost_elements() for PETScWrappers::MPI::BlockVector and TrilinosWrappers::MPI::BlockVector.
    (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; iget_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 +#include + +#include //needed for assembling the matrix using quadrature on each cell +#include +#include +#include +#include +#include +#include +#include //direct solver +#include //definition of the sparse matrix +#include //for the intermediate sparsity pattern structure +#include //CG solver +#include //and a preconditioner +#include //conform hanging nodes DoF to certain constrains to make the solution continuous +#include //need for Schur matrix +#include //triangulation class +#include //standard functions to generate grid +#include +#include +#include //these two to loop over cells/faces +#include // ^ +#include //to use not straightline boundaries +#include //associate DoF to cells/vertices/lines +#include //provides information about the degrees of freedom local to a cell +#include //needed for the creation of sparsity patterns of sparse matrices +#include //use to renumber DoF to have a better sparsity pattern +#include //needed for assembling the matrix using quadrature on each cell +#include +#include +#include +#include +#include +#include //for evaluating solution at a given point +#include /*for vector-valued FE; vector-valude FE will be composed from the regular Q1 included below*/ +#include /*Q1 FE: 1 DoF in each vertex, none on sides...)*/ +#include +#include +#include +#include +#include +#include +#include //sqrt and fabs functions +#include //smart pointers +#include //min,max + + +template +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 dofs_per_block; + + Triangulation triangulation; //a triangulation object of the "dim"-dimensional domain; + hp::DoFHandler dof_handler; //is assotiated with triangulation; + + FESystem elasticity_fe; + FESystem elasticity_w_lagrange_fe; + hp::FECollection 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 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 +SymmetricTensor<4,dim> +get_stress_strain_tensor (const double lambda, const double mu) +{ + SymmetricTensor<4,dim> tmp; + for (unsigned int i=0; i +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function(dim) {} + + virtual void vector_value (const Point &p, + Vector &values) const; + virtual double value (const Point &p, + unsigned int compoment) const; +}; + +template +inline +void BoundaryValues::vector_value (const Point &p, + Vector &values) const +{ + values(0) = -0.001; +} + +template +inline +double BoundaryValues::value (const Point &p, + unsigned int compoment) const +{ + return -0.001; +} + +/// +/* vector valued constrained displacements*/ +/*vector-valued RHS function */ +template +class ConstrainValues : public Function +{ +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 &p, + Vector &values) const;//returns calculated values in the second argument; + virtual void vector_value_list (const std::vector > &points, + std::vector > &value_list) const;//values at several points at once + //prevent from calling virtual function "vector_value" to frequently +}; + +template +ConstrainValues::ConstrainValues () +: +Function (dim) /*pass down to the main class the number of components of the function*/ +{} + +template +inline +void ConstrainValues::vector_value (const Point &p, + Vector &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 +void ConstrainValues::vector_value_list (const std::vector > &points, + std::vector > &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::vector_value (points[p], + value_list[p]); +} + + +/* constructor*/ +template +ElasticProblem::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(degree), dim, //use dim FE_Q of a given degree to represent displacements + FE_Nothing(), dim // zero extension of lagrange multipliers elsewhere in the domain +), +elasticity_w_lagrange_fe( + FE_Q(degree), dim, + FE_Q(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(degree+2)); + quadrature_formula.push_back(QGauss(degree+2)); + + stress_strain_tensor + = get_stress_strain_tensor (/*lambda = */ 9.695e10, + /*mu = */ 7.617e10); +} + +template +ElasticProblem::~ElasticProblem () +{ + dof_handler.clear (); +} + + +template +void ElasticProblem::make_grid () +{ + + Triangulation triangulationL; + Triangulation 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 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;iface(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 +void ElasticProblem::setup_system () +{ + std::vector 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::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 << "; "< 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="< 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 +#include + +#include //triangulation class +#include //standard functions to generate grid +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void test() +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation, -1,0); + triangulation.refine_global (1); + + FESystem elasticity_fe (FE_Q(1), 1, + FE_Nothing(), 1); + FESystem elasticity_w_lagrange_fe(FE_Q(1), 1, + FE_Q(1), 1); + hp::FECollection fe; + fe.push_back (elasticity_fe); + fe.push_back (elasticity_w_lagrange_fe); + + hp::DoFHandler dof_handler (triangulation); + typename hp::DoFHandler::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 -- 2.39.5