Index: include/numerics/vector_value.h =================================================================== --- include/numerics/vector_value.h (revision 5911) +++ include/numerics/vector_value.h (working copy) @@ -181,7 +181,6 @@ { } - #ifdef LIBMESH_USE_COMPLEX_NUMBERS template inline @@ -194,6 +193,12 @@ } #endif +template +inline +Real libmesh_norm(const VectorValue& a) +{return a.size_sq();} + + } // namespace libMesh #endif // #define __vector_value_h__ Index: include/numerics/type_vector.h =================================================================== --- include/numerics/type_vector.h (revision 5911) +++ include/numerics/type_vector.h (working copy) @@ -998,6 +998,11 @@ libmesh_dot(const TypeVector& a, const TypeVector& b) { return a * b; } +template +inline +Real libmesh_norm(const TypeVector& a) +{return a.size_sq();} + } // namespace libMesh #endif // #define __type_vector_h__ Index: include/numerics/raw_accessor.h =================================================================== --- include/numerics/raw_accessor.h (revision 5911) +++ include/numerics/raw_accessor.h (working copy) @@ -54,6 +54,12 @@ typedef Number type; }; + template<> + struct RawFieldType > + { + typedef Number type; + }; + #ifdef LIBMESH_USE_COMPLEX_NUMBERS template <> struct RawFieldType @@ -72,6 +78,12 @@ { typedef Real type; }; + + template<> + struct RawFieldType > + { + typedef Real type; + }; #endif /** @@ -140,7 +152,7 @@ public: typedef TypeNTensor FieldType; - + RawAccessor( FieldType& data, const unsigned int dim ) : _data(data), _dim(dim) @@ -148,16 +160,16 @@ ~RawAccessor(){}; - typename RawFieldType::type& operator()( unsigned int i ) + typename RawFieldType::type& operator()( unsigned int /*i*/ ) { return dummy; } - const typename RawFieldType::type& operator()( unsigned int i ) const + const typename RawFieldType::type& operator()( unsigned int /*i*/ ) const { return dummy; } private: RawAccessor(); - static ScalarType dummy; + ScalarType dummy; FieldType& _data; const unsigned int _dim; Index: include/numerics/function_base.h =================================================================== --- include/numerics/function_base.h (revision 5911) +++ include/numerics/function_base.h (working copy) @@ -48,7 +48,11 @@ * required input of each derived class thwarts the effective * use of the commonly used \p build() member. But afterwards * the virtual members allow the convenient and libMesh-common - * usage through a \p FunctionBase*. + * usage through a \p FunctionBase*. Note that for functor objects + * for vector-valued variables, it is assumed each component is indexed + * contiguously; i.e. if u_var is index 3, then libMesh expects + * the x-component of u_var is index 3, the y-component is index 4, + * and the z-component is index 5. * * @author Daniel Dreyer, 2003 */ Index: include/numerics/type_n_tensor.h =================================================================== --- include/numerics/type_n_tensor.h (revision 5911) +++ include/numerics/type_n_tensor.h (working copy) @@ -152,6 +152,12 @@ contract (const TypeNTensor &) const { return 0; } /** + * Returns the Frobenius norm of the tensor squared, i.e. sum of the + * element magnitudes squared. + */ + Real size_sq() const { return 0.;} + + /** * @returns \p true if two tensors are equal valued. */ bool operator == (const TypeNTensor&) const { return true; } Index: include/error_estimation/exact_solution.h =================================================================== --- include/error_estimation/exact_solution.h (revision 5911) +++ include/error_estimation/exact_solution.h (working copy) @@ -231,6 +231,7 @@ * private function since it is used by the implementation when * solving for several unknowns in several systems. */ + template void _compute_error(const std::string& sys_name, const std::string& unknown_name, std::vector& error_vals); Index: src/error_estimation/exact_solution.C =================================================================== --- src/error_estimation/exact_solution.C (revision 5911) +++ src/error_estimation/exact_solution.C (working copy) @@ -33,6 +33,9 @@ #include "tensor_value.h" #include "vector_value.h" #include "wrapped_function.h" +#include "fe_interface.h" +#include "raw_accessor.h" +#include "type_vector.h" namespace libMesh { @@ -303,9 +306,32 @@ // to the proper location to store the error std::vector& error_vals = this->_check_inputs(sys_name, unknown_name); - this->_compute_error(sys_name, - unknown_name, - error_vals); + + libmesh_assert( _equation_systems.has_system(sys_name) ); + const System& sys = _equation_systems.get_system( sys_name ); + + libmesh_assert( sys.has_variable( unknown_name ) ); + switch( FEInterface::field_type(sys.variable_type( unknown_name )) ) + { + case TYPE_SCALAR: + { + this->_compute_error(sys_name, + unknown_name, + error_vals); + break; + } + case TYPE_VECTOR: + { + this->_compute_error(sys_name, + unknown_name, + error_vals); + break; + } + default: + libmesh_error(); + } + + return; } @@ -453,7 +479,7 @@ - +template< typename OutputShape> void ExactSolution::_compute_error(const std::string& sys_name, const std::string& unknown_name, std::vector& error_vals) @@ -517,19 +543,33 @@ const MeshBase& _mesh = computed_system.get_mesh(); + const unsigned int dim = _mesh.mesh_dimension(); + // Zero the error before summation error_vals = std::vector(5, 0.); // Construct Quadrature rule based on default quadrature order const FEType& fe_type = computed_dof_map.variable_type(var); + unsigned int n_vec_dim = FEInterface::n_vec_dim( _mesh, fe_type ); + + // FIXME: MeshFunction needs to be updated to support vector-valued + // elements before we can use a reference solution. + if( (n_vec_dim > 1) && _equation_systems_fine ) + { + libMesh::err << "Error calculation using reference solution not yet\n" + << "supported for vector-valued elements." + << std::endl; + libmesh_not_implemented(); + } + AutoPtr qrule = - fe_type.default_quadrature_rule (_mesh.mesh_dimension(), + fe_type.default_quadrature_rule (dim, _extra_order); // Construct finite element object - AutoPtr fe(FEBase::build(_mesh.mesh_dimension(), fe_type)); + AutoPtr > fe(FEGenericBase::build(dim, fe_type)); // Attach quadrature rule to FE object fe->attach_quadrature_rule (qrule.get()); @@ -539,14 +579,16 @@ // The value of the shape functions at the quadrature points // i.e. phi(i) = phi_values[i][qp] - const std::vector >& phi_values = fe->get_phi(); + const std::vector >& phi_values = fe->get_phi(); // The value of the shape function gradients at the quadrature points - const std::vector >& dphi_values = fe->get_dphi(); + const std::vector::OutputGradient> >& + dphi_values = fe->get_dphi(); #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES // The value of the shape function second derivatives at the quadrature points - const std::vector >& d2phi_values = fe->get_d2phi(); + const std::vector::OutputTensor> >& + d2phi_values = fe->get_d2phi(); #endif // The XYZ locations (in physical space) of the quadrature points @@ -593,11 +635,11 @@ // Real u_h = 0.; // RealGradient grad_u_h; - Number u_h = 0.; + typename FEGenericBase::OutputNumber u_h = 0.; - Gradient grad_u_h; + typename FEGenericBase::OutputNumberGradient grad_u_h; #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES - Tensor grad2_u_h; + typename FEGenericBase::OutputNumberTensor grad2_u_h; #endif @@ -616,50 +658,82 @@ } // Compute the value of the error at this quadrature point - Number exact_val = 0.0; + typename FEGenericBase::OutputNumber exact_val; + RawAccessor::OutputNumber> exact_val_accessor( exact_val, dim ); if (_exact_values.size() > sys_num && _exact_values[sys_num]) - exact_val = - _exact_values[sys_num]-> - component(var_component, q_point[qp], time); + { + for( unsigned int c = 0; c < n_vec_dim; c++) + exact_val_accessor(c) = + _exact_values[sys_num]-> + component(var_component+c, q_point[qp], time); + } else if (_equation_systems_fine) - exact_val = (*coarse_values)(q_point[qp]); + { + // FIXME: Needs to be updated for vector-valued elements + exact_val = (*coarse_values)(q_point[qp]); + } + const typename FEGenericBase::OutputNumber val_error = u_h - exact_val; - const Number val_error = u_h - exact_val; + // Add the squares of the error to each contribution + Real error_sq = libmesh_norm(val_error); + error_vals[0] += JxW[qp]*error_sq; - // Add the squares of the error to each contribution - error_vals[0] += JxW[qp]*libmesh_norm(val_error); - Real norm = sqrt(libmesh_norm(val_error)); + Real norm = sqrt(error_sq); error_vals[3] += JxW[qp]*norm; + if(error_vals[4]::OutputNumberGradient exact_grad; + RawAccessor::OutputNumberGradient> exact_grad_accessor( exact_grad, dim ); if (_exact_derivs.size() > sys_num && _exact_derivs[sys_num]) - exact_grad = - _exact_derivs[sys_num]-> - component(var_component, q_point[qp], time); + { + for( unsigned int c = 0; c < n_vec_dim; c++) + for( unsigned int d = 0; d < dim; d++ ) + exact_grad_accessor(d + c*dim ) = + _exact_derivs[sys_num]-> + component(var_component+c, q_point[qp], time)(d); + } else if (_equation_systems_fine) - exact_grad = coarse_values->gradient(q_point[qp]); + { + // FIXME: Needs to be updated for vector-valued elements + exact_grad = coarse_values->gradient(q_point[qp]); + } - const Gradient grad_error = grad_u_h - exact_grad; - + const typename FEGenericBase::OutputNumberGradient grad_error = grad_u_h - exact_grad; + error_vals[1] += JxW[qp]*grad_error.size_sq(); - #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES // Compute the value of the error in the hessian at this // quadrature point - Tensor exact_hess; + typename FEGenericBase::OutputNumberTensor exact_hess; + RawAccessor::OutputNumberTensor> exact_hess_accessor( exact_hess, dim ); if (_exact_hessians.size() > sys_num && _exact_hessians[sys_num]) - exact_hess = - _exact_hessians[sys_num]-> - component(var_component, q_point[qp], time); + { + //FIXME: This needs to be implemented to support rank 3 tensors + // which can't happen until type_n_tensor is fully implemented + // and a RawAccessor is fully implemented + if( FEInterface::field_type(fe_type) == TYPE_VECTOR ) + libmesh_not_implemented(); + + for( unsigned int c = 0; c < n_vec_dim; c++) + for( unsigned int d = 0; d < dim; d++ ) + for( unsigned int e =0; e < dim; e++ ) + exact_hess_accessor(d + e*dim + c*dim*dim) = + _exact_hessians[sys_num]-> + component(var_component+c, q_point[qp], time)(d,e); + } else if (_equation_systems_fine) - exact_hess = coarse_values->hessian(q_point[qp]); - - const Tensor grad2_error = grad2_u_h - exact_hess; - + { + // FIXME: Needs to be updated for vector-valued elements + exact_hess = coarse_values->hessian(q_point[qp]); + } + + const typename FEGenericBase::OutputNumberTensor grad2_error = grad2_u_h - exact_hess; + + // FIXME: PB: Is this what we want for rank 3 tensors? error_vals[2] += JxW[qp]*grad2_error.size_sq(); #endif @@ -674,4 +748,8 @@ error_vals[4] = l_infty_norm; } + // Explicit instantiations of templated member functions + template void ExactSolution::_compute_error(const std::string&, const std::string&, std::vector&); + template void ExactSolution::_compute_error(const std::string&, const std::string&, std::vector&); + } // namespace libMesh Index: examples/vector_fe/vector_fe_ex2/laplace_system.C =================================================================== --- examples/vector_fe/vector_fe_ex2/laplace_system.C (revision 5911) +++ examples/vector_fe/vector_fe_ex2/laplace_system.C (working copy) @@ -71,6 +71,7 @@ std::vector vars; vars.push_back( u_var ); + // Note that for vector-valued variables, it is assumed each component is stored contiguously SolutionFunction func( u_var ); this->get_dof_map().add_dirichlet_boundary( libMesh::DirichletBoundary( boundary_ids, vars, &func ) ); Index: examples/vector_fe/vector_fe_ex2/solution_function.h =================================================================== --- examples/vector_fe/vector_fe_ex2/solution_function.h (revision 5911) +++ examples/vector_fe/vector_fe_ex2/solution_function.h (working copy) @@ -40,6 +40,8 @@ { output.zero(); const Real x=p(0), y=p(1), z=p(2); + // libMesh assumes each component of the vector-valued variable is stored + // contiguously. output(_u_var) = soln( 0, x, y, z ); output(_u_var+1) = soln( 1, x, y, z ); output(_u_var+2) = soln( 2, x, y, z ); Index: examples/vector_fe/vector_fe_ex2/vector_fe_ex2.C =================================================================== --- examples/vector_fe/vector_fe_ex2/vector_fe_ex2.C (revision 5911) +++ examples/vector_fe/vector_fe_ex2/vector_fe_ex2.C (working copy) @@ -28,6 +28,7 @@ #include "exodusII_io.h" #include "mesh.h" #include "mesh_generation.h" +#include "exact_solution.h" // The systems and solvers we may use #include "laplace_system.h" @@ -110,11 +111,37 @@ system.solve(); + ExactSolution exact_sol( equation_systems ); + + std::vector* > sols; + std::vector* > grads; + + sols.push_back( new SolutionFunction(system.variable_number("u")) ); + grads.push_back( new SolutionGradient(system.variable_number("u")) ); + + exact_sol.attach_exact_values(sols); + exact_sol.attach_exact_derivs(grads); + + // Use higher quadrature order for more accurate error results + int extra_error_quadrature = infile("extra_error_quadrature",2); + exact_sol.extra_quadrature_order(extra_error_quadrature); + + // Compute the error. + exact_sol.compute_error("Laplace", "u"); + + // Print out the error values + std::cout << "L2-Error is: " + << exact_sol.l2_error("Laplace", "u") + << std::endl; + std::cout << "H1-Error is: " + << exact_sol.h1_error("Laplace", "u") + << std::endl; + #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.