From d34609dd63ec3e4cb7a2b125696221f8ee7fa3c4 Mon Sep 17 00:00:00 2001 From: nivesh Date: Fri, 8 Jun 2018 16:41:52 +0200 Subject: [PATCH] added a test to check correct constraints are used --- tests/fe/fe_enriched_color_08.cc | 285 +++++++++++++++++++++++++++ tests/fe/fe_enriched_color_08.output | 27 +++ 2 files changed, 312 insertions(+) create mode 100644 tests/fe/fe_enriched_color_08.cc create mode 100644 tests/fe/fe_enriched_color_08.output diff --git a/tests/fe/fe_enriched_color_08.cc b/tests/fe/fe_enriched_color_08.cc new file mode 100644 index 0000000000..ec12286211 --- /dev/null +++ b/tests/fe/fe_enriched_color_08.cc @@ -0,0 +1,285 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 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 if correct constraints are assigned for hp::DoFHandler using + * hp::FECollection constructed by ColorEnriched::Helper. + * + * Due to the bug #1496 (https://github.com/dealii/dealii/issues/1496), + * the constraints at certain interfaces may not be correct. + * This tests if special treatment of the interface between enriched cells + * results in the correct constraints despite the bug. + */ + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +#undef DEBUG + +#ifdef DEBUG +// used only for debugging +template +void +plot_shape_function(hp::DoFHandler &dof_handler, unsigned int patches = 5) +{ + deallog << "...start plotting shape function" << std::endl; + deallog << "Patches for output: " << patches << std::endl; + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + // find set of dofs which belong to enriched cells + std::set enriched_cell_dofs; + for (auto cell : dof_handler.active_cell_iterators()) + if (cell->active_fe_index() != 0) + { + unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + std::vector local_dof_indices(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + enriched_cell_dofs.insert(local_dof_indices.begin(), + local_dof_indices.end()); + } + + // output to check if all is good: + std::vector> shape_functions; + std::vector names; + for (auto dof : enriched_cell_dofs) + { + Vector shape_function; + shape_function.reinit(dof_handler.n_dofs()); + shape_function[dof] = 1.0; + + // if the dof is constrained, first output unconstrained vector + names.push_back(std::string("C_") + + dealii::Utilities::int_to_string(dof, 2)); + shape_functions.push_back(shape_function); + + // names.push_back(std::string("UC_") + + // dealii::Utilities::int_to_string(s,2)); + + // // make continuous/constraint: + // constraints.distribute(shape_function); + // shape_functions.push_back(shape_function); + } + + if (dof_handler.n_dofs() < 100) + { + deallog << "...start printing support points" << std::endl; + + std::map> support_points; + MappingQ1 mapping; + hp::MappingCollection hp_mapping; + for (unsigned int i = 0; i < dof_handler.get_fe_collection().size(); ++i) + hp_mapping.push_back(mapping); + DoFTools::map_dofs_to_support_points( + hp_mapping, dof_handler, support_points); + + const std::string base_filename = + "DOFs" + dealii::Utilities::int_to_string(dim) + "_p" + + dealii::Utilities::int_to_string(0); + + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics " + " " + << std::endl + << "unset ytics" << std::endl + << "unset grid" << std::endl + << "unset border" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 " + "offset 1,1 notitle" + << std::endl; + GridOut grid_out; + grid_out.write_gnuplot(dof_handler.get_triangulation(), f); + f << "e" << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, support_points); + + f << "e" << std::endl; + + deallog << "...finished printing support points" << std::endl; + } + + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + + // get material ids: + Vector fe_index(dof_handler.get_triangulation().n_active_cells()); + for (auto cell : dof_handler.active_cell_iterators()) + { + fe_index[cell->active_cell_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(patches); + + std::string filename = "shape_functions.vtu"; + std::ofstream output(filename.c_str()); + data_out.write_vtu(output); + + deallog << "...finished plotting shape functions" << std::endl; +} +#endif + + + +/* + * Testing helper class + */ +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() <= radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + + +/* + * Type used to defined vector of predicates needed by the function + * ColorEnriched::internal::color_predicates. + */ +template +using predicate_function = + std::function::cell_iterator &)>; + + + +int +main(int argc, char **argv) +{ + // Intialize MPI as required by Zoltan library used for graph coloring by this + // test. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + MPILogInitAll all; + + // Make basic grid + const unsigned int dim = 2; + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + Point p1(0, 0), p2(3, 1); + GridGenerator::subdivided_hyper_rectangle(triangulation, {3, 1}, p1, p2); + + // Make predicates resulting in three adjacent domains + // The relevant enrichment functions for the domains are [][1,3][2,3] + // [] means no enrichment function. + std::vector> vec_predicates; + vec_predicates.push_back(EnrichmentPredicate(Point(1.5, 0.5), 0.5)); + vec_predicates.push_back(EnrichmentPredicate(Point(2.5, 0.5), 0.5)); + vec_predicates.push_back(EnrichmentPredicate(Point(2, 0.5), 1)); + + // Construct vector of enrichment functions + std::vector>> vec_enrichments; + vec_enrichments.reserve(vec_predicates.size()); + for (unsigned int i = 0; i < vec_predicates.size(); ++i) + { + // constant function. + ConstantFunction func(10 + i); // constant function + vec_enrichments.push_back(std::make_shared>(func)); + } + + // Construct helper class to construct fe collection + FE_Q fe_base(2); + FE_Q fe_enriched(1); + + static ColorEnriched::Helper fe_space( + fe_base, fe_enriched, vec_predicates, vec_enrichments); + hp::FECollection fe_collection( + fe_space.build_fe_collection(dof_handler)); + + // check if fe_collection is correctly constructed by function + deallog << "fe_collection[index] mapping:" << std::endl; + for (unsigned int index = 0; index != fe_collection.size(); ++index) + { + deallog << "name:" << fe_collection[index].get_name() << std::endl; + deallog << "n_blocks:" << fe_collection[index].n_blocks() << std::endl; + deallog << "n_comp:" << fe_collection[index].n_components() << std::endl; + deallog << "n_dofs:" << fe_collection[index].n_dofs_per_cell() + << std::endl; + } + + deallog << "Dominating set for 1 and 2: " + << fe_collection.find_least_face_dominating_fe({1, 2}) << std::endl; + + dof_handler.distribute_dofs(fe_collection); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + constraints.print(deallog.get_file_stream()); + +#ifdef DEBUG + plot_shape_function(dof_handler); +#endif + + dof_handler.clear(); + return 0; +} diff --git a/tests/fe/fe_enriched_color_08.output b/tests/fe/fe_enriched_color_08.output new file mode 100644 index 0000000000..675f3567bb --- /dev/null +++ b/tests/fe/fe_enriched_color_08.output @@ -0,0 +1,27 @@ + +DEAL:0::fe_collection[index] mapping: +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)] +DEAL:0::n_blocks:4 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:9 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)-FE_Q<2>(1)] +DEAL:0::n_blocks:4 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:17 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)-FE_Q<2>(1)] +DEAL:0::n_blocks:4 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:17 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)-FE_Q<2>(1)] +DEAL:0::n_blocks:4 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:13 +DEAL:0::Dominating set for 1 and 2: 3 + 9 = 0 + 10 = 0 + 12 = 0 + 14 = 0 + 15 = 0 + 17 = 0 + 23 = 0 + 27 = 0 -- 2.39.5