From fd26695f99680c18cd22a99f325742571d69fd2e Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 1 Sep 2015 20:52:44 +0200 Subject: [PATCH] test make_hanging_node_constraints() with p-ref and neither_element_dominates --- .../hp/hp_constraints_neither_dominate_01.cc | 260 ++++++++++++++++++ .../hp_constraints_neither_dominate_01.output | 101 +++++++ 2 files changed, 361 insertions(+) create mode 100644 tests/hp/hp_constraints_neither_dominate_01.cc create mode 100644 tests/hp/hp_constraints_neither_dominate_01.output diff --git a/tests/hp/hp_constraints_neither_dominate_01.cc b/tests/hp/hp_constraints_neither_dominate_01.cc new file mode 100644 index 0000000000..36824cecce --- /dev/null +++ b/tests/hp/hp_constraints_neither_dominate_01.cc @@ -0,0 +1,260 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 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 hp_constraints when neither_element_dominates. +// 2D and 3D p-refinement only for 2 elements. +// Check continuity across the face by evaluating the +// field from both sides. + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +#include +#include + +// #define DEBUG_OUTPUT_VTK + +const double eps = 1e-10; + +using namespace dealii; + +template +void test2cells(const unsigned int p1=2, + const unsigned int p2=1) +{ + Triangulation triangulation; + { + Triangulation triangulationL; + Triangulation triangulationR; + GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain + GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain + Point shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulationR); + GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation); + } + + hp::DoFHandler dof_handler(triangulation); + + hp::FECollection fe_collection; + fe_collection.push_back(FESystem(FE_Q(p1),1, + FE_Nothing(1,true),1)); + fe_collection.push_back(FESystem(FE_Q(p2),1, + FE_Q(1),1)); + // push back to be able to resolve hp constrains no matter what: + fe_collection.push_back(FESystem(FE_Q(p2),1, + FE_Nothing(1,true),1)); + + hp::QCollection< dim-1 > q_face_collection; + q_face_collection.push_back(QIterated(QTrapez<1>(),2)); + q_face_collection.push_back(QIterated(QTrapez<1>(),2)); + q_face_collection.push_back(QIterated(QTrapez<1>(),2)); + + deallog << "2cells: "<set_active_fe_index(1);//Q(p2)xQ(1) + + dof_handler.distribute_dofs(fe_collection); + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + constraints.print(deallog.get_file_stream()); + +#ifdef DEBUG_OUTPUT_VTK + // output to check if all is good: + std::vector> shape_functions; + std::vector names; + for (unsigned int s=0; s < dof_handler.n_dofs(); s++) + { + Vector shape_function; + shape_function.reinit(dof_handler.n_dofs()); + shape_function[s] = 1.0; + + // if the dof is constrained, first output unconstrained vector + if (constraints.is_constrained(s)) + { + names.push_back(std::string("UN_") + + dealii::Utilities::int_to_string(s,2)); + shape_functions.push_back(shape_function); + } + + names.push_back(std::string("N_") + + dealii::Utilities::int_to_string(s,2)); + + // make continuous/constrain: + constraints.distribute(shape_function); + shape_functions.push_back(shape_function); + } + + DataOut> data_out; + data_out.attach_dof_handler (dof_handler); + + // get material ids: + Vector fe_index(triangulation.n_active_cells()); + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (unsigned int index=0; cell!=endc; ++cell,++index) + { + fe_index[index] = cell->active_fe_index(); + } + data_out.add_data_vector(fe_index, "fe_index"); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + data_out.add_data_vector (shape_functions[i], names[i]); + + data_out.build_patches(15); + + std::string filename = "2cell-shape_functions_p1="+std::to_string(p1) + +"_p2="+std::to_string(p2) + +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + // fill some vector + Vector solution(dof_handler.n_dofs()); + for (unsigned int dof=0; dof < dof_handler.n_dofs(); dof++) + solution[dof] = 21.0*(dof%2) + 0.5 + dof%3; + + constraints.distribute(solution); + + + // evaluate field at the interface: + hp::FEFaceValues< dim > fe_face_values_hp(fe_collection, + q_face_collection, + update_values | update_quadrature_points); + + std::vector local_face_dof_indices; + std::vector > values; + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + { + const unsigned int fe_index = cell->active_fe_index(); + local_face_dof_indices.resize(fe_collection[fe_index].dofs_per_face); + for (unsigned int f=0; f::faces_per_cell; ++f) + if (!cell->at_boundary(f)) // that's enough for our simple mesh + { + deallog << "cell="<(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing + test2cells<2>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing + test2cells<3>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing + test2cells<3>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/hp/hp_constraints_neither_dominate_01.output b/tests/hp/hp_constraints_neither_dominate_01.output new file mode 100644 index 0000000000..10f1da2b1b --- /dev/null +++ b/tests/hp/hp_constraints_neither_dominate_01.output @@ -0,0 +1,101 @@ + +DEAL::2cells: FESystem<2>[FE_Q<2>(1)-FE_Nothing<2>()] vs FESystem<2>[FE_Q<2>(2)-FE_Q<2>(1)] + 2 = 0 + 5 = 0 + 7 11: 0.500000 + 7 13: 0.500000 +DEAL::cell=0.0 face=1 +DEAL::u[0.0000 -1.0000]={23.5000,0} +DEAL::u[0.0000 -0.5000]={23.0000,0} +DEAL::u[0.0000 0.0000]={22.5000,0} +DEAL::cell=0.1 face=0 +DEAL::u[0.0000 -1.0000]={23.5000,0} +DEAL::u[0.0000 -0.5000]={23.0000,0} +DEAL::u[0.0000 0.0000]={22.5000,0} +DEAL::2cells: FESystem<2>[FE_Q<2>(2)-FE_Nothing<2>()] vs FESystem<2>[FE_Q<2>(1)-FE_Q<2>(1)] + 2 = 0 + 5 = 0 + 10 6: 0.500000 + 10 8: 0.500000 +DEAL::cell=0.0 face=1 +DEAL::u[0.0000 -1.0000]={0.5000,0} +DEAL::u[0.0000 -0.5000]={1.5000,0} +DEAL::u[0.0000 0.0000]={2.5000,0} +DEAL::cell=0.1 face=0 +DEAL::u[0.0000 -1.0000]={0.5000,0} +DEAL::u[0.0000 -0.5000]={1.5000,0} +DEAL::u[0.0000 0.0000]={2.5000,0} +DEAL::2cells: FESystem<3>[FE_Q<3>(1)-FE_Nothing<3>()] vs FESystem<3>[FE_Q<3>(2)-FE_Q<3>(1)] + 2 = 0 + 5 = 0 + 8 = 0 + 11 = 0 + 13 31: 0.500000 + 13 33: 0.500000 + 17 35: 0.500000 + 17 37: 0.500000 + 21 31: 0.500000 + 21 35: 0.500000 + 23 33: 0.500000 + 23 37: 0.500000 + 25 31: 0.250000 + 25 33: 0.250000 + 25 35: 0.250000 + 25 37: 0.250000 +DEAL::cell=0.0 face=1 +DEAL::u[0.0000 -1.0000 -1.0000]={22.5000,0} +DEAL::u[0.0000 -0.5000 -1.0000]={22.0000,0} +DEAL::u[0.0000 0.0000 -1.0000]={21.5000,0} +DEAL::u[0.0000 -1.0000 -0.5000]={23.0000,0} +DEAL::u[0.0000 -0.5000 -0.5000]={22.5000,0} +DEAL::u[0.0000 0.0000 -0.5000]={22.0000,0} +DEAL::u[0.0000 -1.0000 0.0000]={23.5000,0} +DEAL::u[0.0000 -0.5000 0.0000]={23.0000,0} +DEAL::u[0.0000 0.0000 0.0000]={22.5000,0} +DEAL::cell=0.1 face=0 +DEAL::u[0.0000 -1.0000 -1.0000]={22.5000,0} +DEAL::u[0.0000 -0.5000 -1.0000]={22.0000,0} +DEAL::u[0.0000 0.0000 -1.0000]={21.5000,0} +DEAL::u[0.0000 -1.0000 -0.5000]={23.0000,0} +DEAL::u[0.0000 -0.5000 -0.5000]={22.5000,0} +DEAL::u[0.0000 0.0000 -0.5000]={22.0000,0} +DEAL::u[0.0000 -1.0000 0.0000]={23.5000,0} +DEAL::u[0.0000 -0.5000 0.0000]={23.0000,0} +DEAL::u[0.0000 0.0000 0.0000]={22.5000,0} +DEAL::2cells: FESystem<3>[FE_Q<3>(2)-FE_Nothing<3>()] vs FESystem<3>[FE_Q<3>(1)-FE_Q<3>(1)] + 2 = 0 + 5 = 0 + 8 = 0 + 11 = 0 + 20 12: 0.500000 + 20 14: 0.500000 + 24 16: 0.500000 + 24 18: 0.500000 + 28 12: 0.500000 + 28 16: 0.500000 + 30 14: 0.500000 + 30 18: 0.500000 + 32 12: 0.250000 + 32 14: 0.250000 + 32 16: 0.250000 + 32 18: 0.250000 +DEAL::cell=0.0 face=1 +DEAL::u[0.0000 -1.0000 -1.0000]={0.5000,0} +DEAL::u[0.0000 -0.5000 -1.0000]={1.5000,0} +DEAL::u[0.0000 0.0000 -1.0000]={2.5000,0} +DEAL::u[0.0000 -1.0000 -0.5000]={1.0000,0} +DEAL::u[0.0000 -0.5000 -0.5000]={1.2500,0} +DEAL::u[0.0000 0.0000 -0.5000]={1.5000,0} +DEAL::u[0.0000 -1.0000 0.0000]={1.5000,0} +DEAL::u[0.0000 -0.5000 0.0000]={1.0000,0} +DEAL::u[0.0000 0.0000 0.0000]={0.5000,0} +DEAL::cell=0.1 face=0 +DEAL::u[0.0000 -1.0000 -1.0000]={0.5000,0} +DEAL::u[0.0000 -0.5000 -1.0000]={1.5000,0} +DEAL::u[0.0000 0.0000 -1.0000]={2.5000,0} +DEAL::u[0.0000 -1.0000 -0.5000]={1.0000,0} +DEAL::u[0.0000 -0.5000 -0.5000]={1.2500,0} +DEAL::u[0.0000 0.0000 -0.5000]={1.5000,0} +DEAL::u[0.0000 -1.0000 0.0000]={1.5000,0} +DEAL::u[0.0000 -0.5000 0.0000]={1.0000,0} +DEAL::u[0.0000 0.0000 0.0000]={0.5000,0} -- 2.39.5