diff --git a/examples/fem_system/fem_system_ex4/Makefile.am b/examples/fem_system/fem_system_ex4/Makefile.am new file mode 100644 index 00000000000..1e918b3d2c4 --- /dev/null +++ b/examples/fem_system/fem_system_ex4/Makefile.am @@ -0,0 +1,21 @@ +example_name = fem_system_ex4 +check_SCRIPTS = run.sh +install_dir = $(examples_install_path)/fem_system/ex1 +data = fem_system_ex4.C fem_system_ex4.in heatsystem.C heatsystem.h run.sh +sources = $(data) run.sh + +CLEANFILES = out_*.e + +# also need to link files for VPATH builds +if LIBMESH_VPATH_BUILD + BUILT_SOURCES = .linkstamp +.linkstamp: + -rm -f fem_system_ex4.in && $(LN_S) -f $(srcdir)/fem_system_ex4.in . + $(AM_V_GEN)touch .linkstamp + + CLEANFILES += fem_system_ex4.in .linkstamp +endif + +############################################## +# include common example environment +include $(top_srcdir)/examples/Make.common diff --git a/examples/fem_system/fem_system_ex4/fem_system_ex4.C b/examples/fem_system/fem_system_ex4/fem_system_ex4.C new file mode 100644 index 00000000000..8d1b3e18fa0 --- /dev/null +++ b/examples/fem_system/fem_system_ex4/fem_system_ex4.C @@ -0,0 +1,302 @@ +/* The libMesh Finite Element Library. */ +/* Copyright (C) 2003 Benjamin S. Kirk */ + +/* This library is free software; you can 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. */ + +/* This library is distributed in the hope that it will be useful, */ +/* but WITHOUT ANY WARRANTY; without even the implied warranty of */ +/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */ +/* Lesser General Public License for more details. */ + +/* You should have received a copy of the GNU Lesser General Public */ +/* License along with this library; if not, write to the Free Software */ +/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ + +//

FEMSystem Example 4 - Mixed-dimension heat transfer equation +// with FEMSystem

+// +// This example shows how elements of multiple dimensions can be +// linked and computed upon simultaneously within the +// DifferentiableSystem class framework + +// C++ includes +#include + +// Basic include files +#include "libmesh/equation_systems.h" +#include "libmesh/error_vector.h" +#include "libmesh/exact_solution.h" +#include "libmesh/getpot.h" +#include "libmesh/gmv_io.h" +#include "libmesh/exodusII_io.h" +#include "libmesh/kelly_error_estimator.h" +#include "libmesh/mesh.h" +#include "libmesh/mesh_generation.h" +#include "libmesh/mesh_refinement.h" +#include "libmesh/parsed_function.h" +#include "libmesh/uniform_refinement_estimator.h" + +// The systems and solvers we may use +#include "heatsystem.h" +#include "libmesh/diff_solver.h" +#include "libmesh/euler_solver.h" +#include "libmesh/steady_solver.h" + +// Bring in everything from the libMesh namespace +using namespace libMesh; + +// The main program. +int main (int argc, char** argv) +{ + // Initialize libMesh. + LibMeshInit init (argc, argv); + +#ifndef LIBMESH_ENABLE_AMR + libmesh_example_requires(false, "--enable-amr"); +#else + + // Parse the input file + GetPot infile("fem_system_ex4.in"); + + // Read in parameters from the input file + const Real global_tolerance = infile("global_tolerance", 0.); + const unsigned int nelem_target = infile("n_elements", 400); + const Real deltat = infile("deltat", 0.005); + const unsigned int coarsegridsize = infile("coarsegridsize", 20); + const unsigned int coarserefinements = infile("coarserefinements", 0); + const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10); + const unsigned int dim = infile("dimension", 2); + + // Skip higher-dimensional examples on a lower-dimensional libMesh build + libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support"); + + // We have only defined 2 and 3 dimensional problems + libmesh_assert (dim == 2 || dim == 3); + + // Create a mesh, with dimension to be overridden later, distributed + // across the default MPI communicator. + Mesh mesh(init.comm()); + + // And an object to refine it + MeshRefinement mesh_refinement(mesh); + mesh_refinement.coarsen_by_parents() = true; + mesh_refinement.absolute_global_tolerance() = global_tolerance; + mesh_refinement.nelem_target() = nelem_target; + mesh_refinement.refine_fraction() = 0.3; + mesh_refinement.coarsen_fraction() = 0.3; + mesh_refinement.coarsen_threshold() = 0.1; + + // Use the MeshTools::Generation mesh generator to create a uniform + // grid on the square or cube. We crop the domain at y=2/3 to allow + // for a homogeneous Neumann BC in our benchmark there. + boundary_id_type bcid = 3; // +y in 3D + if (dim == 2) + { + MeshTools::Generation::build_square + (mesh, + coarsegridsize, + coarsegridsize*2/3, // int arithmetic best we can do here + 0., 1., + 0., 2./3., + QUAD9); + bcid = 2; // +y in 2D + } + else if (dim == 3) + { + MeshTools::Generation::build_cube + (mesh, + coarsegridsize, + coarsegridsize*2/3, + coarsegridsize, + 0., 1., + 0., 2./3., + 0., 1., + HEX27); + } + + { + // Add boundary elements corresponding to the +y boundary of our + // volume mesh + std::set bcids; + bcids.insert(bcid); + mesh.get_boundary_info().add_elements(bcids, mesh); + mesh.prepare_for_use(); + } + + mesh_refinement.uniformly_refine(coarserefinements); + + // Print information about the mesh to the screen. + mesh.print_info(); + + // Create an equation systems object. + EquationSystems equation_systems (mesh); + + // Declare the system "Heat" and its variables. + HeatSystem & system = + equation_systems.add_system ("Heat"); + + // Solve this as a steady system + system.time_solver = + UniquePtr(new SteadySolver(system)); + + // Initialize the system + equation_systems.init (); + + // Set the time stepping options + system.deltat = deltat; + + // And the nonlinear solver options + DiffSolver &solver = *(system.time_solver->diff_solver().get()); + solver.quiet = infile("solver_quiet", true); + solver.verbose = !solver.quiet; + solver.max_nonlinear_iterations = + infile("max_nonlinear_iterations", 15); + solver.relative_step_tolerance = + infile("relative_step_tolerance", 1.e-3); + solver.relative_residual_tolerance = + infile("relative_residual_tolerance", 0.0); + solver.absolute_residual_tolerance = + infile("absolute_residual_tolerance", 0.0); + + // And the linear solver options + solver.max_linear_iterations = + infile("max_linear_iterations", 50000); + solver.initial_linear_tolerance = + infile("initial_linear_tolerance", 1.e-3); + + // Print information about the system to the screen. + equation_systems.print_info(); + + // Adaptively solve the steady solution + unsigned int a_step = 0; + for (; a_step != max_adaptivesteps; ++a_step) + { + system.solve(); + + system.postprocess(); + + ErrorVector error; + + UniquePtr error_estimator; + + // To solve to a tolerance in this problem we + // need a better estimator than Kelly + if (global_tolerance != 0.) + { + // We can't adapt to both a tolerance and a mesh + // size at once + libmesh_assert_equal_to (nelem_target, 0); + + UniformRefinementEstimator *u = + new UniformRefinementEstimator; + + // The lid-driven cavity problem isn't in H1, so + // lets estimate L2 error + u->error_norm = L2; + + error_estimator.reset(u); + } + else + { + // If we aren't adapting to a tolerance we need a + // target mesh size + libmesh_assert_greater (nelem_target, 0); + + // Kelly is a lousy estimator to use for a problem + // not in H1 - if we were doing more than a few + // timesteps we'd need to turn off or limit the + // maximum level of our adaptivity eventually + error_estimator.reset(new KellyErrorEstimator); + } + + error_estimator->estimate_error(system, error); + + // Print out status at each adaptive step. + Real global_error = error.l2_norm(); + std::cout << "Adaptive step " << a_step << ": " << std::endl; + if (global_tolerance != 0.) + std::cout << "Global_error = " << global_error + << std::endl; + if (global_tolerance != 0.) + std::cout << "Worst element error = " << error.maximum() + << ", mean = " << error.mean() << std::endl; + + if (global_tolerance != 0.) + { + // If we've reached our desired tolerance, we + // don't need any more adaptive steps + if (global_error < global_tolerance) + break; + mesh_refinement.flag_elements_by_error_tolerance(error); + } + else + { + // If flag_elements_by_nelem_target returns true, this + // should be our last adaptive step. + if (mesh_refinement.flag_elements_by_nelem_target(error)) + { + mesh_refinement.refine_and_coarsen_elements(); + equation_systems.reinit(); + a_step = max_adaptivesteps; + break; + } + } + + // Carry out the adaptive mesh refinement/coarsening + mesh_refinement.refine_and_coarsen_elements(); + equation_systems.reinit(); + + std::cout << "Refined mesh to " + << mesh.n_active_elem() + << " active elements and " + << equation_systems.n_active_dofs() + << " active dofs." << std::endl; + } + // Do one last solve if necessary + if (a_step == max_adaptivesteps) + { + system.solve(); + + system.postprocess(); + } + + +#ifdef LIBMESH_HAVE_EXODUS_API + // We want to write the file in the ExodusII format, but we don't + // yet support mixed-dimension meshes there? +/* ExodusII_IO(mesh).write_equation_systems + ("out.e", equation_systems); +*/ +#endif // #ifdef LIBMESH_HAVE_EXODUS_API + +#ifdef LIBMESH_HAVE_GMV + GMVIO(mesh).write_equation_systems + ("out.gmv", equation_systems); +#endif // #ifdef LIBMESH_HAVE_GMV + +#ifdef LIBMESH_HAVE_FPARSER + // Check that we got close to the analytic solution + ExactSolution exact_sol(equation_systems); + const std::string exact_str = (dim == 2) ? + "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)"; + ParsedFunction exact_func(exact_str); + exact_sol.attach_exact_value(0, &exact_func); + exact_sol.compute_error("Heat", "T"); + + Number err = exact_sol.l2_error("Heat", "T"); + + // Print out the error value + std::cout << "L2-Error is: " << err << std::endl; + + libmesh_assert_less(err, 2e-3); + +#endif // #ifdef LIBMESH_HAVE_FPARSER + +#endif // #ifndef LIBMESH_ENABLE_AMR + + // All done. + return 0; +} diff --git a/examples/fem_system/fem_system_ex4/fem_system_ex4.in b/examples/fem_system/fem_system_ex4/fem_system_ex4.in new file mode 100644 index 00000000000..f9eaa61aa5d --- /dev/null +++ b/examples/fem_system/fem_system_ex4/fem_system_ex4.in @@ -0,0 +1,26 @@ +# The global FEM error tolerance at each timestep +# Make this nonzero to solve to a specified tolerance +# This will probably break with KellyErrorIndicator +# const Real global_tolerance = 1.e-3; +global_tolerance = 0.0 + +# The desired number of active mesh elements +# Make this nonzero to solve to a specified mesh size +n_elements = 400 + +# The coarse grid size from which to start adaptivity +coarsegridsize = 20 + +# The maximum number of adaptive steps per timestep +max_adaptivesteps = 0 + +# Turn this off to see the NewtonSolver chatter +solver_quiet = false + +# Solve the 2D or 3D problem +dimension = 2 + +# Nonlinear solver tolerance +relative_step_tolerance = 1e-9 + +relative_residual_tolerance = 1.e-9 diff --git a/examples/fem_system/fem_system_ex4/heatsystem.C b/examples/fem_system/fem_system_ex4/heatsystem.C new file mode 100644 index 00000000000..2965801220b --- /dev/null +++ b/examples/fem_system/fem_system_ex4/heatsystem.C @@ -0,0 +1,143 @@ +#include "heatsystem.h" + +#include "libmesh/dirichlet_boundaries.h" +#include "libmesh/dof_map.h" +#include "libmesh/fe_base.h" +#include "libmesh/fe_interface.h" +#include "libmesh/fem_context.h" +#include "libmesh/getpot.h" +#include "libmesh/mesh.h" +#include "libmesh/quadrature.h" +#include "libmesh/string_to_enum.h" +#include "libmesh/zero_function.h" + +using namespace libMesh; + +void HeatSystem::init_data () +{ + T_var = this->add_variable + ("T", static_cast(_fe_order), + Utility::string_to_enum(_fe_family)); + + const unsigned int dim = this->get_mesh().mesh_dimension(); + + // Add dirichlet boundaries on all but the boundary element side + const boundary_id_type all_ids[6] = {0, 1, 2, 3, 4, 5}; + std::set nonyplus_bdys(all_ids, all_ids+(dim*2)); + const boundary_id_type yplus_id = (dim == 3) ? 3 : 2; + nonyplus_bdys.erase(yplus_id); + + std::vector T_only(1, T_var); + ZeroFunction zero; + + this->get_dof_map().add_dirichlet_boundary + (DirichletBoundary (nonyplus_bdys, T_only, &zero)); + + // Do the parent's initialization after variables are defined + FEMSystem::init_data(); + + this->time_evolving(0); +} + + + +void HeatSystem::init_context(DiffContext &context) +{ + FEMContext &c = libmesh_cast_ref(context); + + const std::set& elem_dims = + c.elem_dimensions(); + + for (std::set::const_iterator dim_it = + elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it) + { + const unsigned char dim = *dim_it; + + FEBase* fe = NULL; + + c.get_element_fe( T_var, fe, dim ); + + fe->get_JxW(); // For integration + fe->get_dphi(); // For bilinear form + fe->get_xyz(); // For forcing + fe->get_phi(); // For forcing + } + + FEMSystem::init_context(context); +} + + +bool HeatSystem::element_time_derivative (bool request_jacobian, + DiffContext &context) +{ + FEMContext &c = libmesh_cast_ref(context); + + const unsigned int mesh_dim = + c.get_system().get_mesh().mesh_dimension(); + + // First we get some references to cell-specific data that + // will be used to assemble the linear system. + const unsigned int dim = c.get_elem().dim(); + FEBase* fe = NULL; + c.get_element_fe(T_var, fe, dim); + + // Element Jacobian * quadrature weights for interior integration + const std::vector &JxW = fe->get_JxW(); + + const std::vector &xyz = fe->get_xyz(); + + const std::vector > &phi = fe->get_phi(); + + const std::vector > &dphi = fe->get_dphi(); + + // The number of local degrees of freedom in each variable + const unsigned int n_T_dofs = c.get_dof_indices(T_var).size(); + + // The subvectors and submatrices we need to fill: + DenseSubMatrix &K = c.get_elem_jacobian(T_var, T_var); + DenseSubVector &F = c.get_elem_residual(T_var); + + // Now we will build the element Jacobian and residual. + // Constructing the residual requires the solution and its + // gradient from the previous timestep. This must be + // calculated at each quadrature point by summing the + // solution degree-of-freedom values by the appropriate + // weight functions. + unsigned int n_qpoints = c.get_element_qrule().n_points(); + + for (unsigned int qp=0; qp != n_qpoints; qp++) + { + // Compute the solution gradient at the Newton iterate + Gradient grad_T = c.interior_gradient(T_var, qp); + + const Number k = _k[dim]; + + const Point &p = xyz[qp]; + + // solution + laplacian depend on problem dimension + const Number u_exact = (mesh_dim == 2) ? + std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) : + std::sin(libMesh::pi*p(0)) * std::sin(libMesh::pi*p(1)) * + std::sin(libMesh::pi*p(2)); + + // Only apply forcing to interior elements + const Number forcing = (dim == mesh_dim) ? + -k * dim * libMesh::pi * libMesh::pi * u_exact : 0; + + const Number JxWxNK = JxW[qp] * -k; + + for (unsigned int i=0; i != n_T_dofs; i++) + F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]); + if (request_jacobian) + { + const Number JxWxNKxD = JxWxNK * + context.get_elem_solution_derivative(); + + for (unsigned int i=0; i != n_T_dofs; i++) + for (unsigned int j=0; j != n_T_dofs; ++j) + K(i,j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]); + } + } // end of the quadrature point qp-loop + + return request_jacobian; +} diff --git a/examples/fem_system/fem_system_ex4/heatsystem.h b/examples/fem_system/fem_system_ex4/heatsystem.h new file mode 100644 index 00000000000..a8b53f10c99 --- /dev/null +++ b/examples/fem_system/fem_system_ex4/heatsystem.h @@ -0,0 +1,48 @@ +#include "libmesh/enum_fe_family.h" +#include "libmesh/fem_system.h" + +// FEMSystem, TimeSolver and NewtonSolver will handle most tasks, +// but we must specify element residuals +class HeatSystem : public libMesh::FEMSystem +{ +public: + // Constructor + HeatSystem(libMesh::EquationSystems& es, + const std::string& name, + const unsigned int number) + : libMesh::FEMSystem(es, name, number), + _fe_family("LAGRANGE"), _fe_order(1) { + // Get the conductivity ratios right for both 2D and 3D + // benchmarks + _k[1] = 1/libMesh::pi/std::sqrt(libMesh::Real(3)); + _k[2] = 1; + _k[3] = 2*libMesh::pi*std::sqrt(libMesh::Real(3)); + + } + + std::string & fe_family() { return _fe_family; } + unsigned int & fe_order() { return _fe_order; } + +protected: + // System initialization + virtual void init_data (); + + // Context initialization + virtual void init_context (libMesh::DiffContext &context); + + // Element residual and jacobian calculations + // Time dependent parts + virtual bool element_time_derivative (bool request_jacobian, + libMesh::DiffContext &context); + + // The conductivity for the various dimensional elements, indexed by + // dim (with _k[0] unused) for simplicity + libMesh::Real _k[4]; + + // The FE type to use + std::string _fe_family; + unsigned int _fe_order; + + // The variable index (yes, this will be 0...) + unsigned int T_var; +}; diff --git a/examples/fem_system/fem_system_ex4/run.sh b/examples/fem_system/fem_system_ex4/run.sh new file mode 100755 index 00000000000..00402f41dfe --- /dev/null +++ b/examples/fem_system/fem_system_ex4/run.sh @@ -0,0 +1,11 @@ +#!/bin/bash + +#set -x + +source $LIBMESH_DIR/examples/run_common.sh + +example_name=fem_system_ex4 + +example_dir=examples/fem_system/$example_name + +run_example "$example_name"