Index: examples/vector_fe/vector_fe_ex2/laplace_system.C =================================================================== --- examples/vector_fe/vector_fe_ex2/laplace_system.C (revision 0) +++ examples/vector_fe/vector_fe_ex2/laplace_system.C (working copy) @@ -0,0 +1,251 @@ +/* The Next Great 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 */ + +#include "getpot.h" + +#include "laplace_system.h" + +#include "boundary_info.h" +#include "dof_map.h" +#include "fe_base.h" +#include "fe_interface.h" +#include "fem_context.h" +#include "mesh.h" +#include "quadrature.h" +#include "string_to_enum.h" + + +// Bring in everything from the libMesh namespace +using namespace libMesh; + +LaplaceSystem::LaplaceSystem( EquationSystems& es, + const std::string& name, + const unsigned int number) + : FEMSystem(es, name, number) +{ + return; +} + +void LaplaceSystem::init_data () +{ + // Check the input file for Reynolds number, application type, + // approximation type + GetPot infile("laplace.in"); + + // Add the solution variable + u_var = this->add_variable ("u", SECOND, LAGRANGE_VEC); + + this->time_evolving(u_var); + + // Useful debugging options + // Set verify_analytic_jacobians to 1e-6 to use + this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.); + this->print_jacobians = infile("print_jacobians", false); + this->print_element_jacobians = infile("print_element_jacobians", false); + + // Do the parent's initialization after variables and boundary constraints are defined + FEMSystem::init_data(); +} + + + +void LaplaceSystem::init_context(DiffContext &context) +{ + FEMContext &c = libmesh_cast_ref(context); + + // Get finite element object + FEGenericBase* fe; + c.get_element_fe( u_var, fe ); + + // We should prerequest all the data + // we will need to build the linear system. + fe->get_JxW(); + fe->get_phi(); + fe->get_dphi(); + fe->get_xyz(); + + FEGenericBase* side_fe; + c.get_side_fe( u_var, side_fe ); + + side_fe->get_JxW(); + side_fe->get_phi(); +} + + +bool LaplaceSystem::element_time_derivative (bool request_jacobian, + DiffContext &context) +{ + FEMContext &c = libmesh_cast_ref(context); + + // Get finite element object + FEGenericBase* fe = NULL; + c.get_element_fe( u_var, fe ); + + // First we get some references to cell-specific data that + // will be used to assemble the linear system. + + // Element Jacobian * quadrature weights for interior integration + const std::vector &JxW = fe->get_JxW(); + + // The velocity shape functions at interior quadrature points. + const std::vector >& phi = fe->get_phi(); + + // The velocity shape function gradients at interior + // quadrature points. + const std::vector >& grad_phi = fe->get_dphi(); + + const std::vector& qpoint = fe->get_xyz(); + + // The number of local degrees of freedom in each variable + const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); + + DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; + + DenseSubVector &Fu = *c.elem_subresiduals[u_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. + const unsigned int n_qpoints = c.element_qrule->n_points(); + + for (unsigned int qp=0; qp != n_qpoints; qp++) + { + RealTensor grad_u; + + c.interior_gradient( u_var, qp, grad_u ); + + + // Value of the forcing function at this quadrature point + RealGradient f = this->forcing(qpoint[qp]); + + // First, an i-loop over the velocity degrees of freedom. + // We know that n_u_dofs == n_v_dofs so we can compute contributions + // for both at the same time. + for (unsigned int i=0; i != n_u_dofs; i++) + { + Fu(i) += ( grad_u.contract(grad_phi[i][qp]) - f*phi[i][qp] )*JxW[qp]; + + if (request_jacobian) + { + // Matrix contributions for the uu and vv couplings. + for (unsigned int j=0; j != n_u_dofs; j++) + { + Kuu(i,j) += grad_phi[j][qp].contract(grad_phi[i][qp])*JxW[qp]; + + } + } + } + } // end of the quadrature point qp-loop + + return request_jacobian; +} + +bool LaplaceSystem::side_constraint (bool request_jacobian, + DiffContext &context) +{ + FEMContext &c = libmesh_cast_ref(context); + + // Get finite element object + FEGenericBase* side_fe = NULL; + c.get_side_fe( u_var, side_fe ); + + // First we get some references to cell-specific data that + // will be used to assemble the linear system. + + // Element Jacobian * quadrature weights for interior integration + const std::vector &JxW = side_fe->get_JxW(); + + // The velocity shape functions at interior quadrature points. + const std::vector >& phi = side_fe->get_phi(); + + // The number of local degrees of freedom in each variable + const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); + + const std::vector& qpoint = side_fe->get_xyz(); + + // The penalty value. \frac{1}{\epsilon} + // in the discussion above. + const Real penalty = 1.e10; + + DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; + DenseSubVector &Fu = *c.elem_subresiduals[u_var]; + + const unsigned int n_qpoints = c.side_qrule->n_points(); + + for (unsigned int qp=0; qp != n_qpoints; qp++) + { + RealGradient u; + c.side_value( u_var, qp, u ); + + RealGradient u_exact( this->exact_solution( 0, qpoint[qp](0), qpoint[qp](1) ), + this->exact_solution( 1, qpoint[qp](0), qpoint[qp](1) )); + + for (unsigned int i=0; i != n_u_dofs; i++) + { + Fu(i) += penalty*(u - u_exact)*phi[i][qp]*JxW[qp]; + + if (request_jacobian) + { + for (unsigned int j=0; j != n_u_dofs; j++) + Kuu(i,j) += penalty*phi[j][qp]*phi[i][qp]*JxW[qp]; + } + } + } + + return request_jacobian; +} + + +Real LaplaceSystem::exact_solution(unsigned int component, + Real x, Real y, Real z) +{ + switch(component) + { + case 0: + return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z); + + case 1: + return sin(.5*pi*x)*cos(.5*pi*y)*cos(.5*pi*z); + + default: + libmesh_error(); + } +} + +RealGradient LaplaceSystem::forcing( const Point& p ) +{ + Real x = p(0); Real y = p(1); + + const Real eps = 1.e-3; + + const Real fx = -(exact_solution(0,x,y-eps) + + exact_solution(0,x,y+eps) + + exact_solution(0,x-eps,y) + + exact_solution(0,x+eps,y) - + 4.*exact_solution(0,x,y))/eps/eps; + + const Real fy = -(exact_solution(1,x,y-eps) + + exact_solution(1,x,y+eps) + + exact_solution(1,x-eps,y) + + exact_solution(1,x+eps,y) - + 4.*exact_solution(1,x,y))/eps/eps; + + return RealGradient( fx, fy ); +} Index: examples/vector_fe/vector_fe_ex2/laplace.in =================================================================== --- examples/vector_fe/vector_fe_ex2/laplace.in (revision 0) +++ examples/vector_fe/vector_fe_ex2/laplace.in (working copy) @@ -0,0 +1,5 @@ + +verify_analytic_jacobians = 0. +print_jacobians = false +print_element_jacobians = false + Index: examples/vector_fe/vector_fe_ex2/vector_fe_ex2.C =================================================================== --- examples/vector_fe/vector_fe_ex2/vector_fe_ex2.C (revision 0) +++ examples/vector_fe/vector_fe_ex2/vector_fe_ex2.C (working copy) @@ -0,0 +1,120 @@ +/* The Next Great 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 1 - Unsteady Navier-Stokes Equations with + // FEMSystem

+ // + // This example shows how the transient nonlinear problem from + // example 13 can be solved using the + // DifferentiableSystem class framework + +// Basic include files +#include "equation_systems.h" +#include "getpot.h" +#include "exodusII_io.h" +#include "mesh.h" +#include "mesh_generation.h" + +// The systems and solvers we may use +#include "laplace_system.h" +#include "diff_solver.h" +#include "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); + + // Parse the input file + GetPot infile("vector_fe_ex2.in"); + + // Read in parameters from the input file + const unsigned int grid_size = infile( "grid_size", 2 ); + + // Skip higher-dimensional examples on a lower-dimensional libMesh build + libmesh_example_assert(2 <= LIBMESH_DIM, "2D/3D support"); + + // Create a mesh. + Mesh mesh; + + // Use the MeshTools::Generation mesh generator to create a uniform + // grid on the square [-1,1]^D. We instruct the mesh generator + // to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 + // elements in 3D. Building these higher-order elements allows + // us to use higher-order approximation, as in example 3. + MeshTools::Generation::build_square (mesh, + grid_size, + grid_size, + -1., 1., + -1., 1., + QUAD9); + + // Print information about the mesh to the screen. + mesh.print_info(); + + // Create an equation systems object. + EquationSystems equation_systems (mesh); + + // Declare the system "Navier-Stokes" and its variables. + LaplaceSystem & system = + equation_systems.add_system ("Laplace"); + + // This example only implements the steady-state problem + system.time_solver = + AutoPtr(new SteadySolver(system)); + + // Initialize the system + equation_systems.init (); + + // 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(); + + system.solve(); + +#ifdef LIBMESH_HAVE_EXODUS_API + + // We write the file in the ExodusII format. + ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems); + +#endif // #ifdef LIBMESH_HAVE_EXODUS_API + + // All done. + return 0; +} Index: examples/vector_fe/vector_fe_ex2/laplace_system.h =================================================================== --- examples/vector_fe/vector_fe_ex2/laplace_system.h (revision 0) +++ examples/vector_fe/vector_fe_ex2/laplace_system.h (working copy) @@ -0,0 +1,61 @@ +/* The Next Great 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 */ + +// DiffSystem framework files +#include "fem_system.h" +#include "vector_value.h" +#include "tensor_value.h" + +using namespace libMesh; + +// FEMSystem, TimeSolver and NewtonSolver will handle most tasks, +// but we must specify element residuals +class LaplaceSystem : public FEMSystem +{ + +public: + +// Constructor +LaplaceSystem( EquationSystems& es, + const std::string& name, + const unsigned int number); + +// System initialization +virtual void init_data (); + +// Context initialization +virtual void init_context(DiffContext &context); + +// Element residual and jacobian calculations +// Time dependent parts +virtual bool element_time_derivative (bool request_jacobian, + DiffContext& context); + +// Constraint part +virtual bool side_constraint (bool request_jacobian, + DiffContext& context); + +protected: +// Indices for each variable; +unsigned int u_var; + +// Returns the value of a forcing function at point p. +RealGradient forcing(const Point& p); + +Real exact_solution( unsigned int c, + Real x, Real y, Real z = 0.0 ); +}; Index: examples/vector_fe/vector_fe_ex2/Makefile =================================================================== --- examples/vector_fe/vector_fe_ex2/Makefile (revision 0) +++ examples/vector_fe/vector_fe_ex2/Makefile (working copy) @@ -0,0 +1,74 @@ +# The location of the mesh library +LIBMESH_DIR ?= ../../.. + +# include the library options determined by configure. This will +# set the variables INCLUDE and LIBS that we will need to build and +# link with the library. +include $(LIBMESH_DIR)/Make.common + + +############################################################################### +# File management. This is where the source, header, and object files are +# defined + +# +# source files +srcfiles := $(wildcard *.C) + +# +# object files +objects := $(patsubst %.C, %.$(obj-suffix), $(srcfiles)) +############################################################################### + + + +.PHONY: clean clobber distclean + +############################################################################### +# Target: +# +target := ./vector_fe_ex2-$(METHOD) + + +all:: $(target) + +# Production rules: how to make the target - depends on library configuration +$(target): $(objects) + @echo "Linking "$@"..." + @$(libmesh_CXX) $(libmesh_CPPFLAGS) $(libmesh_CXXFLAGS) $(objects) -o $@ $(libmesh_LIBS) $(libmesh_LDFLAGS) + +# Useful rules. +clean: + @rm -f $(objects) *~ .depend + +clobber: + @$(MAKE) clean + @rm -f $(target) out.e + +distclean: + @$(MAKE) clobber + @rm -f *.o *.g.o *.pg.o .depend + +run: $(target) + @echo "***************************************************************" + @echo "* Running Example " $(LIBMESH_RUN) $(target) $(LIBMESH_OPTIONS) + @echo "***************************************************************" + @echo " " + @$(LIBMESH_RUN) $(target) $(LIBMESH_OPTIONS) + @echo " " + @echo "***************************************************************" + @echo "* Done Running Example " $(LIBMESH_RUN) $(target) $(LIBMESH_OPTIONS) + @echo "***************************************************************" + + +# include the dependency list +include .depend + + +# +# Dependencies +# +.depend: $(srcfiles) $(LIBMESH_DIR)/include/*/*.h + @$(perl) $(LIBMESH_DIR)/contrib/bin/make_dependencies.pl -I. $(foreach i, $(wildcard $(LIBMESH_DIR)/include/*), -I$(i)) "-S\$$(obj-suffix)" $(srcfiles) > .depend + +############################################################################### Index: examples/vector_fe/vector_fe_ex2/vector_fe_ex2.in =================================================================== --- examples/vector_fe/vector_fe_ex2/vector_fe_ex2.in (revision 0) +++ examples/vector_fe/vector_fe_ex2/vector_fe_ex2.in (working copy) @@ -0,0 +1,11 @@ +# The coarse grid size from which to start adaptivity +grid_size = 20 + +# 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-3 Index: src/systems/fem_context.C =================================================================== --- src/systems/fem_context.C (revision 5729) +++ src/systems/fem_context.C (working copy) @@ -32,10 +32,6 @@ namespace libMesh { - - - - FEMContext::FEMContext (const System &sys) : DiffContext(sys), element_qrule(NULL), side_qrule(NULL), @@ -77,31 +73,63 @@ (1, sys.extra_quadrature_order).release(); // Next, create finite element objects + // Preserving backward compatibility here for now + // Should move to the protected/FEAbstract interface element_fe_var.resize(n_vars); side_fe_var.resize(n_vars); if (dim == 3) edge_fe_var.resize(n_vars); + _element_fe_var.resize(n_vars); + _side_fe_var.resize(n_vars); + if (dim == 3) + _edge_fe_var.resize(n_vars); + for (unsigned int i=0; i != n_vars; ++i) { FEType fe_type = sys.variable_type(i); - if (element_fe[fe_type] == NULL) - { - element_fe[fe_type] = FEBase::build(dim, fe_type).release(); - element_fe[fe_type]->attach_quadrature_rule(element_qrule); - side_fe[fe_type] = FEBase::build(dim, fe_type).release(); - side_fe[fe_type]->attach_quadrature_rule(side_qrule); - if (dim == 3) - { - edge_fe[fe_type] = FEBase::build(dim, fe_type).release(); - edge_fe[fe_type]->attach_quadrature_rule(edge_qrule); - } - } - element_fe_var[i] = element_fe[fe_type]; - side_fe_var[i] = side_fe[fe_type]; + // Preserving backward compatibility here for now + // Should move to the protected/FEAbstract interface + if( FEInterface::field_type( fe_type ) == TYPE_SCALAR ) + { + if( element_fe[fe_type] == NULL ) + { + element_fe[fe_type] = FEBase::build(dim, fe_type).release(); + element_fe[fe_type]->attach_quadrature_rule(element_qrule); + side_fe[fe_type] = FEBase::build(dim, fe_type).release(); + side_fe[fe_type]->attach_quadrature_rule(side_qrule); + + if (dim == 3) + { + edge_fe[fe_type] = FEBase::build(dim, fe_type).release(); + edge_fe[fe_type]->attach_quadrature_rule(edge_qrule); + } + } + element_fe_var[i] = element_fe[fe_type]; + side_fe_var[i] = side_fe[fe_type]; + if (dim == 3) + edge_fe_var[i] = edge_fe[fe_type]; + } + + if ( _element_fe[fe_type] == NULL ) + { + _element_fe[fe_type] = FEAbstract::build(dim, fe_type).release(); + _element_fe[fe_type]->attach_quadrature_rule(element_qrule); + _side_fe[fe_type] = FEAbstract::build(dim, fe_type).release(); + _side_fe[fe_type]->attach_quadrature_rule(side_qrule); + + if (dim == 3) + { + _edge_fe[fe_type] = FEAbstract::build(dim, fe_type).release(); + _edge_fe[fe_type]->attach_quadrature_rule(edge_qrule); + } + } + _element_fe_var[i] = _element_fe[fe_type]; + _side_fe_var[i] = _side_fe[fe_type]; if (dim == 3) - edge_fe_var[i] = edge_fe[fe_type]; + _edge_fe_var[i] = _edge_fe[fe_type]; + } } @@ -110,6 +138,8 @@ { // We don't want to store AutoPtrs in STL containers, but we don't // want to leak memory either + // Preserving backward compatibility here for now + // Should move to the protected/FEAbstract interface for (std::map::iterator i = element_fe.begin(); i != element_fe.end(); ++i) delete i->second; @@ -125,6 +155,23 @@ delete i->second; edge_fe.clear(); + + for (std::map::iterator i = _element_fe.begin(); + i != _element_fe.end(); ++i) + delete i->second; + _element_fe.clear(); + + for (std::map::iterator i = _side_fe.begin(); + i != _side_fe.end(); ++i) + delete i->second; + _side_fe.clear(); + + for (std::map::iterator i = _edge_fe.begin(); + i != _edge_fe.end(); ++i) + delete i->second; + _edge_fe.clear(); + + delete element_qrule; element_qrule = NULL; @@ -162,8 +209,36 @@ return u; } +template +void FEMContext::interior_value(unsigned int var, unsigned int qp, OutputShape& u) const +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + // Get finite element object + FEGenericBase* fe = NULL; + this->get_element_fe( var, fe ); + + // Get shape function values at quadrature point + const std::vector > &phi = fe->get_phi(); + + // Accumulate solution value + u = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + u += phi[l][qp] * coef(l); + + return; +} + + + Gradient FEMContext::interior_gradient(unsigned int var, unsigned int qp) const { // Get local-to-global dof index lookup @@ -188,8 +263,37 @@ return du; } +template +void FEMContext::interior_gradient(unsigned int var, unsigned int qp, + typename FEGenericBase::OutputGradient & du) const +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + // Get finite element object + FEGenericBase* fe = NULL; + this->get_element_fe( var, fe ); + + // Get shape function values at quadrature point + const std::vector::OutputGradient> > &dphi = fe->get_dphi(); + + // Accumulate solution derivatives + du = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + du.add_scaled(dphi[l][qp], coef(l)); + + return; +} + + + #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES Tensor FEMContext::interior_hessian(unsigned int var, unsigned int qp) const { @@ -243,7 +347,36 @@ } +template +void FEMContext::side_value(unsigned int var, unsigned int qp, OutputShape& u) const +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get finite element object + FEGenericBase* side_fe = NULL; + this->get_side_fe( var, side_fe ); + + // Get shape function values at quadrature point + const std::vector > &phi = side_fe->get_phi(); + + // Accumulate solution value + u = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + u += phi[l][qp] * coef(l); + + return; +} + + + Gradient FEMContext::side_gradient(unsigned int var, unsigned int qp) const { // Get local-to-global dof index lookup @@ -269,7 +402,36 @@ } +template +void FEMContext::side_gradient(unsigned int var, unsigned int qp, OutputGradient& du) const +{ + // Get local-to-global dof index lookup + libmesh_assert (dof_indices.size() > var); + const unsigned int n_dofs = dof_indices_var[var].size(); + // Get current local coefficients + libmesh_assert (elem_subsolutions.size() > var); + libmesh_assert (elem_subsolutions[var] != NULL); + DenseSubVector &coef = *elem_subsolutions[var]; + + // Get finite element object + FEGenericBase* side_fe = NULL; + this->get_side_fe( var, side_fe ); + + // Get shape function values at quadrature point + const std::vector > &dphi = side_fe->get_dphi(); + + // Accumulate solution derivatives + du = 0.; + + for (unsigned int l=0; l != n_dofs; l++) + du.add_scaled(dphi[l][qp], coef(l)); + + return; +} + + + #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES Tensor FEMContext::side_hessian(unsigned int var, unsigned int qp) const { @@ -729,6 +891,14 @@ { i->second->reinit(elem); } + + + std::map::iterator local_fe_end = _element_fe.end(); + for (std::map::iterator i = _element_fe.begin(); + i != local_fe_end; ++i) + { + i->second->reinit(elem); + } } @@ -742,6 +912,13 @@ { i->second->reinit(elem, side); } + + std::map::iterator local_fe_end = _side_fe.end(); + for (std::map::iterator i = _side_fe.begin(); + i != local_fe_end; ++i) + { + i->second->reinit(elem, side); + } } @@ -758,6 +935,13 @@ { i->second->edge_reinit(elem, edge); } + + std::map::iterator local_fe_end = _edge_fe.end(); + for (std::map::iterator i = _edge_fe.begin(); + i != local_fe_end; ++i) + { + i->second->edge_reinit(elem, edge); + } } @@ -971,5 +1155,17 @@ this->time = theta*(this->system_time + deltat) + (1.-theta)*this->system_time; } +// Instantiate member function templates +template void FEMContext::interior_value(unsigned int, unsigned int, Number&) const; +template void FEMContext::interior_value(unsigned int, unsigned int, Gradient&) const; + template void FEMContext::interior_gradient(unsigned int, unsigned int, FEGenericBase::OutputGradient&) const; +template void FEMContext::interior_gradient(unsigned int, unsigned int, FEGenericBase::OutputGradient&) const; + +template void FEMContext::side_value(unsigned int, unsigned int, Number&) const; +template void FEMContext::side_value(unsigned int, unsigned int, Gradient&) const; + +template void FEMContext::side_gradient(unsigned int, unsigned int, FEGenericBase::OutputGradient&) const; +template void FEMContext::side_gradient(unsigned int, unsigned int, FEGenericBase::OutputGradient&) const; + } // namespace libMesh Index: include/systems/fem_context.h =================================================================== --- include/systems/fem_context.h (revision 5729) +++ include/systems/fem_context.h (working copy) @@ -1,4 +1,3 @@ - // The libMesh Finite Element Library. // Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner @@ -25,6 +24,7 @@ #include "diff_context.h" #include "vector_value.h" #include "fe_type.h" +#include "fe_base.h" #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES #include "tensor_value.h" @@ -186,6 +186,27 @@ */ Tensor fixed_point_hessian (unsigned int var, const Point &p) const; + template + void get_element_fe( unsigned int var, FEGenericBase *& fe ) const; + + template + void get_side_fe( unsigned int var, FEGenericBase *& fe ) const; + + template + void get_edge_fe( unsigned int var, FEGenericBase *& fe ) const; + + template + void interior_value(unsigned int var, unsigned int qp, OutputShape& u) const; + + template + void interior_gradient(unsigned int var, unsigned int qp, typename FEGenericBase::OutputGradient & du) const; + + template + void side_value(unsigned int var, unsigned int qp, OutputShape& u) const; + + template + void side_gradient(unsigned int var, unsigned int qp, OutputGradient& du) const; + #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES /** @@ -316,6 +337,24 @@ */ unsigned char dim; +protected: + + /** + * Finite element objects for each variable's interior, sides and edges. + */ + std::map _element_fe; + std::map _side_fe; + std::map _edge_fe; + + + /** + * Pointers to the same finite element objects, but indexed + * by variable number + */ + std::vector _element_fe_var; + std::vector _side_fe_var; + std::vector _edge_fe_var; + private: /** * Uses the coordinate data specified by mesh_*_position configuration @@ -340,8 +379,6 @@ // ------------------------------------------------------------ // FEMContext inline methods - - inline void FEMContext::elem_position_set(Real theta) { @@ -349,6 +386,31 @@ this->_do_elem_position_set(theta); } +template +inline +void FEMContext::get_element_fe( unsigned int var, FEGenericBase *& fe ) const +{ + libmesh_assert( var < _element_fe_var.size() ); + fe = libmesh_cast_ptr*>( _element_fe_var[var] ); +} + +template +inline +void FEMContext::get_side_fe( unsigned int var, FEGenericBase *& fe ) const +{ + libmesh_assert( var < _side_fe_var.size() ); + fe = libmesh_cast_ptr*>( _side_fe_var[var] ); +} + +template +inline +void FEMContext::get_edge_fe( unsigned int var, FEGenericBase *& fe ) const +{ + libmesh_assert( var < _edge_fe_var.size() ); + fe = libmesh_cast_ptr*>( _edge_fe_var[var] ); +} + + } // namespace libMesh -#endif +#endif //__fem_context_h__