From aeab9cd1f2d6a88cb7702a60093a0f5aaea35daa Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 12 Dec 2014 13:05:30 -0600 Subject: [PATCH] Add another testcase. This is for the same problem as reported by Krishna Garikipati on the mailing list this week. He verified that this used to fail and works now, so I would like to add it to the testsuite to make sure we don't regress in this functionality. --- tests/hp/interpolate_nothing_02.cc | 230 +++++++++++++++++++++++++ tests/hp/interpolate_nothing_02.output | 4 + 2 files changed, 234 insertions(+) create mode 100644 tests/hp/interpolate_nothing_02.cc create mode 100644 tests/hp/interpolate_nothing_02.output diff --git a/tests/hp/interpolate_nothing_02.cc b/tests/hp/interpolate_nothing_02.cc new file mode 100644 index 0000000000..bcacb1105e --- /dev/null +++ b/tests/hp/interpolate_nothing_02.cc @@ -0,0 +1,230 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 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. +// +// --------------------------------------------------------------------- + + + +// like _01, but with the real testcase by Krishna Garikipati + +#include "../tests.h" +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +#define DIMS 3 +#define totalDOF 4 + +template +class diffusionMechanics +{ +public: + diffusionMechanics (const unsigned int mech_degree, const unsigned int diff_degree); + ~diffusionMechanics(); + void run (); + + //Input lengths of rectangular prism + double alen, blen, clen; +private: + enum + {omega1_domain_id, + omega2_domain_id + }; + + static bool + cell_is_in_omega1_domain (const typename hp::DoFHandler::cell_iterator &cell); + static bool + cell_is_in_omega2_domain (const typename hp::DoFHandler::cell_iterator &cell); + void set_active_fe_indices (); + void setup_dofs (); + void setup_system(); + + const unsigned int mech_degree; + const unsigned int diff_degree; + + Triangulation triangulation; + FESystem omega1_fe; + FESystem omega2_fe; + hp::FECollection fe_collection; + hp::DoFHandler dof_handler; + QGauss quadrature_formula; + QGauss face_quadrature_formula; + SparsityPattern sparsity_pattern; + std::map boundary_values; + + Vector system_rhs, U, Un, dU, U0; + + //solution variables + std::vector nodal_solution_names; +}; + +// Constructor +template +diffusionMechanics::diffusionMechanics (const unsigned int mech_degree, const unsigned int diff_degree): +mech_degree (mech_degree), +diff_degree (diff_degree), +omega1_fe (FE_Q(mech_degree), dim, + FE_Q(diff_degree), 1), +omega2_fe (FE_Q(mech_degree), dim, + FE_Nothing(), 1), +dof_handler (triangulation), quadrature_formula(3), face_quadrature_formula(2){ + + //Nodal Solution names + for (unsigned int i=0; i +diffusionMechanics::~diffusionMechanics (){dof_handler.clear ();} + +//Initial conditions - use only if FE_Nothing not used. Else an error occurs with VectorTools::interpolate +template +class InitialConditions: public Function{ +public: + InitialConditions (): Function(totalDOF){} + void vector_value (const Point &p, Vector &values) const{ + Assert (values.size() == totalDOF, ExcDimensionMismatch (values.size(), totalDOF)); + values(totalDOF-4)=0; // u=0 + values(totalDOF-3)=0; + values(totalDOF-2)=0; + values(totalDOF-1)=1.0; + }; +}; +// Boolean functions for specifying sub-domains +template +bool +diffusionMechanics:: +cell_is_in_omega1_domain (const typename hp::DoFHandler::cell_iterator &cell) +{ + return(cell->material_id() == omega1_domain_id); +} +template +bool +diffusionMechanics:: +cell_is_in_omega2_domain (const typename hp::DoFHandler::cell_iterator &cell) +{ + return(cell->material_id() == omega2_domain_id); +} + +// Set active fe indices in each sub-domain +template +void +diffusionMechanics::set_active_fe_indices() +{ + for (typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + { + if (cell_is_in_omega1_domain(cell)) + cell->set_active_fe_index(0); + else if (cell_is_in_omega2_domain(cell)) + cell->set_active_fe_index(1); + else + Assert(false, ExcNotImplemented()); + } +} + +// Call function to set active indices and distribute dofs +template +void +diffusionMechanics::setup_dofs() +{ + set_active_fe_indices (); + dof_handler.distribute_dofs(fe_collection); +} + +//====================done==================================== + +//Setup +template +void diffusionMechanics::setup_system(){ + dof_handler.distribute_dofs (fe_collection); + sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); + U.reinit (dof_handler.n_dofs()); Un.reinit (dof_handler.n_dofs()); dU.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); U0.reinit (dof_handler.n_dofs()); + deallog << " Number of active cells: " << triangulation.n_active_cells() << std::endl; + deallog << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; +} + + +//Run +template +void diffusionMechanics::run (){ + + //subdivided_hyper + std::vector repetitions(3); + repetitions[0] = 2; + repetitions[1] = 2; + repetitions[2] = 4; + GridGenerator::subdivided_hyper_rectangle (triangulation, repetitions, + Point<3>(0.0,0.0,0.0), + Point<3>(alen,blen,clen)); + + //Mark cells by sub-domain + for (typename Triangulation::active_cell_iterator cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell->center()[dim-1] >= 0.75*clen) + cell->set_material_id(omega1_domain_id); + else + cell->set_material_id(omega2_domain_id); + setup_dofs(); + setup_system(); + + //Initial conditions + VectorTools::interpolate(dof_handler, InitialConditions(), Un); + // VectorTools::interpolate(dof_handler, ZeroFunction(fe_collection.n_components()), Un); + U=Un; U0=Un; + + deallog << Un.l2_norm() << std::endl; +} + + +int main () +{ + std::ofstream logfile("output"); + logfile.precision (3); + deallog << std::setprecision(3); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + diffusionMechanics problem(1,1); + problem.run (); +} + diff --git a/tests/hp/interpolate_nothing_02.output b/tests/hp/interpolate_nothing_02.output new file mode 100644 index 0000000000..d8064e07ae --- /dev/null +++ b/tests/hp/interpolate_nothing_02.output @@ -0,0 +1,4 @@ + +DEAL:: Number of active cells: 16 +DEAL:: Number of degrees of freedom: 153 +DEAL::4.24 -- 2.39.5