From 11c9f16691c70dd912550f46c2c524d01e41fe7b Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 28 Aug 2012 07:38:02 +0000 Subject: [PATCH] failing constraints git-svn-id: https://svn.dealii.org/trunk@26138 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/fail/constraints_c1.cc | 157 +++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 tests/fail/constraints_c1.cc diff --git a/tests/fail/constraints_c1.cc b/tests/fail/constraints_c1.cc new file mode 100644 index 0000000000..dcec5fdb84 --- /dev/null +++ b/tests/fail/constraints_c1.cc @@ -0,0 +1,157 @@ +/********************************************************************** + * $Id: constraint_test.cc -1 $ + * + * Copyright Guido Kanschat, 2010, 2012 + * + **********************************************************************/ + +// H2-elements implemented through constraints on the degrees of +// freedom. After adding all constrained lines, the program hangs in +// constraints.close(). + + +#include +#include + +#include +#include + +#include +#include +#include + +using namespace dealii; + +template +void +setup_constraints(const DoFHandler& dof_handler) +{ + ConstraintMatrix constraints; + constraints.clear(); + const FiniteElement& fe = dof_handler.get_fe(); + + // Set up derivative constraints to + // make element C1 + std::vector > > vertex_constraints( + GeometryInfo::vertices_per_cell, std::vector >( + dim+1, std::vector(fe.dofs_per_cell, 0.))); + + for (unsigned int vertex=0;vertex::vertices_per_cell;++vertex) + { + Point v; + for (unsigned int d=0;d 1.e-12) + vertex_constraints[vertex][d][j] = f; + } + const double f = fe.shape_grad_grad(j,v)[0][1]; + if (std::fabs(f) > 1.e-12) + vertex_constraints[vertex][dim][j] = f; + } + } + + std::vector cell_indices(fe.dofs_per_cell); + std::vector neighbor_indices(fe.dofs_per_cell); + + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); + cell != dof_handler.end();++cell) + { + deallog << "New cell" << std::endl; + + cell->get_dof_indices(cell_indices); + + // Do lower left and upper + // right vertex + for (unsigned int face=0;face::faces_per_cell;face+=2) + { + if (cell->at_boundary(face)) continue; + for (unsigned int fvertex=0;fvertex::vertices_per_face;++fvertex) + { + const unsigned int other_face = (face==0) ? (2+fvertex) : fvertex; + if (cell->at_boundary(other_face)) continue; + + const unsigned int vertex = GeometryInfo::face_to_cell_vertices(face, fvertex); + typename DoFHandler::active_cell_iterator neighbor = cell->neighbor(other_face); + const unsigned int neighbor_face = cell->neighbor_of_neighbor(other_face); + neighbor->get_dof_indices(neighbor_indices); + unsigned int neighbor_vertex = GeometryInfo::face_to_cell_vertices(neighbor_face, fvertex); + + const unsigned int d = GeometryInfo::unit_normal_direction[other_face]; + + std::vector > rhs; + const unsigned int constrained = fe.face_to_cell_index( + (fvertex == 0) ? 2 : fe.dofs_per_face-1, face); + double constrained_weight = 0.; + + for (unsigned int i=0;i +void +run(const FiniteElement& fe) +{ + deallog << "Element: " << fe.get_name() << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global(2); + deallog << "Triangulation " + << triangulation.n_active_cells() << " cells, " + << triangulation.n_levels() << " levels" << std::endl; + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + setup_constraints (dof_handler); +} + + +int main() +{ + const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output"); + std::ofstream logfile(logname.c_str()); + deallog.attach(logfile); + + FE_Q<2> fe(3); + run(fe); +} -- 2.39.5