Index: examples/fem_system/fem_system_ex1/naviersystem.C =================================================================== --- examples/fem_system/fem_system_ex1/naviersystem.C (revision 6229) +++ examples/fem_system/fem_system_ex1/naviersystem.C (working copy) @@ -179,18 +179,26 @@ { FEMContext &c = libmesh_cast_ref(context); + FEBase* u_elem_fe; + FEBase* p_elem_fe; + FEBase* u_side_fe; + + c.get_element_fe( u_var, u_elem_fe ); + c.get_element_fe( p_var, p_elem_fe ); + c.get_side_fe( u_var, u_side_fe ); + // We should prerequest all the data // we will need to build the linear system. - c.element_fe_var[u_var]->get_JxW(); - c.element_fe_var[u_var]->get_phi(); - c.element_fe_var[u_var]->get_dphi(); - c.element_fe_var[u_var]->get_xyz(); + u_elem_fe->get_JxW(); + u_elem_fe->get_phi(); + u_elem_fe->get_dphi(); + u_elem_fe->get_xyz(); - c.element_fe_var[p_var]->get_phi(); + p_elem_fe->get_phi(); - c.side_fe_var[u_var]->get_JxW(); - c.side_fe_var[u_var]->get_phi(); - c.side_fe_var[u_var]->get_xyz(); + u_side_fe->get_JxW(); + u_side_fe->get_phi(); + u_side_fe->get_xyz(); } @@ -199,53 +207,54 @@ { FEMContext &c = libmesh_cast_ref(context); + FEBase* u_elem_fe; + FEBase* p_elem_fe; + + c.get_element_fe( u_var, u_elem_fe ); + c.get_element_fe( p_var, p_elem_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 = - c.element_fe_var[u_var]->get_JxW(); + const std::vector &JxW = u_elem_fe->get_JxW(); // The velocity shape functions at interior quadrature points. - const std::vector >& phi = - c.element_fe_var[u_var]->get_phi(); + const std::vector >& phi = u_elem_fe->get_phi(); // The velocity shape function gradients at interior // quadrature points. - const std::vector >& dphi = - c.element_fe_var[u_var]->get_dphi(); + const std::vector >& dphi = u_elem_fe->get_dphi(); // The pressure shape functions at interior // quadrature points. - const std::vector >& psi = - c.element_fe_var[p_var]->get_phi(); + const std::vector >& psi = p_elem_fe->get_phi(); // Physical location of the quadrature points - const std::vector& qpoint = - c.element_fe_var[u_var]->get_xyz(); - + const std::vector& qpoint = u_elem_fe->get_xyz(); + // The number of local degrees of freedom in each variable - const unsigned int n_p_dofs = c.dof_indices_var[p_var].size(); - const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); - libmesh_assert_equal_to (n_u_dofs, c.dof_indices_var[v_var].size()); + const unsigned int n_p_dofs = c.get_dof_indices( p_var ).size(); + const unsigned int n_u_dofs = c.get_dof_indices( u_var ).size(); + libmesh_assert_equal_to (n_u_dofs, c.get_dof_indices( v_var ).size()); // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); - DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; - DenseSubMatrix &Kvv = *c.elem_subjacobians[v_var][v_var]; - DenseSubMatrix &Kww = *c.elem_subjacobians[w_var][w_var]; - DenseSubMatrix &Kuv = *c.elem_subjacobians[u_var][v_var]; - DenseSubMatrix &Kuw = *c.elem_subjacobians[u_var][w_var]; - DenseSubMatrix &Kvu = *c.elem_subjacobians[v_var][u_var]; - DenseSubMatrix &Kvw = *c.elem_subjacobians[v_var][w_var]; - DenseSubMatrix &Kwu = *c.elem_subjacobians[w_var][u_var]; - DenseSubMatrix &Kwv = *c.elem_subjacobians[w_var][v_var]; - DenseSubMatrix &Kup = *c.elem_subjacobians[u_var][p_var]; - DenseSubMatrix &Kvp = *c.elem_subjacobians[v_var][p_var]; - DenseSubMatrix &Kwp = *c.elem_subjacobians[w_var][p_var]; - DenseSubVector &Fu = *c.elem_subresiduals[u_var]; - DenseSubVector &Fv = *c.elem_subresiduals[v_var]; - DenseSubVector &Fw = *c.elem_subresiduals[w_var]; + DenseSubMatrix &Kuu = c.get_elem_jacobian( u_var, u_var ); + DenseSubMatrix &Kvv = c.get_elem_jacobian( v_var, v_var ); + DenseSubMatrix &Kww = c.get_elem_jacobian( w_var, w_var ); + DenseSubMatrix &Kuv = c.get_elem_jacobian( u_var, v_var ); + DenseSubMatrix &Kuw = c.get_elem_jacobian( u_var, w_var ); + DenseSubMatrix &Kvu = c.get_elem_jacobian( v_var, u_var ); + DenseSubMatrix &Kvw = c.get_elem_jacobian( v_var, w_var ); + DenseSubMatrix &Kwu = c.get_elem_jacobian( w_var, u_var ); + DenseSubMatrix &Kwv = c.get_elem_jacobian( w_var, v_var ); + DenseSubMatrix &Kup = c.get_elem_jacobian( u_var, p_var ); + DenseSubMatrix &Kvp = c.get_elem_jacobian( v_var, p_var ); + DenseSubMatrix &Kwp = c.get_elem_jacobian( w_var, p_var ); + DenseSubVector &Fu = c.get_elem_residual( u_var ); + DenseSubVector &Fv = c.get_elem_residual( v_var ); + DenseSubVector &Fw = c.get_elem_residual( w_var ); // Now we will build the element Jacobian and residual. // Constructing the residual requires the solution and its @@ -378,32 +387,36 @@ { FEMContext &c = libmesh_cast_ref(context); + FEBase* u_elem_fe; + FEBase* p_elem_fe; + + c.get_element_fe( u_var, u_elem_fe ); + c.get_element_fe( p_var, p_elem_fe ); + // Here we define some references to cell-specific data that // will be used to assemble the linear system. // Element Jacobian * quadrature weight for interior integration - const std::vector &JxW = c.element_fe_var[u_var]->get_JxW(); + const std::vector &JxW = u_elem_fe->get_JxW(); // The velocity shape function gradients at interior // quadrature points. - const std::vector >& dphi = - c.element_fe_var[u_var]->get_dphi(); + const std::vector >& dphi = u_elem_fe->get_dphi(); // The pressure shape functions at interior // quadrature points. - const std::vector >& psi = - c.element_fe_var[p_var]->get_phi(); + const std::vector >& psi = p_elem_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 unsigned int n_p_dofs = c.dof_indices_var[p_var].size(); + const unsigned int n_u_dofs = c.get_dof_indices( u_var ).size(); + const unsigned int n_p_dofs = c.get_dof_indices( p_var ).size(); // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); - DenseSubMatrix &Kpu = *c.elem_subjacobians[p_var][u_var]; - DenseSubMatrix &Kpv = *c.elem_subjacobians[p_var][v_var]; - DenseSubMatrix &Kpw = *c.elem_subjacobians[p_var][w_var]; - DenseSubVector &Fp = *c.elem_subresiduals[p_var]; + DenseSubMatrix &Kpu = c.get_elem_jacobian( p_var, u_var ); + DenseSubMatrix &Kpv = c.get_elem_jacobian( p_var, v_var ); + DenseSubMatrix &Kpw = c.get_elem_jacobian( p_var, w_var ); + DenseSubVector &Fp = c.get_elem_residual( p_var ); // Add the constraint given by the continuity equation unsigned int n_qpoints = (c.get_element_qrule())->n_points(); @@ -424,9 +437,9 @@ Fp(i) += JxW[qp] * psi[i][qp] * (grad_w(2)); - if (request_jacobian && c.elem_solution_derivative) + if (request_jacobian && c.get_elem_solution_derivative()) { - libmesh_assert_equal_to (c.elem_solution_derivative, 1.0); + libmesh_assert_equal_to (c.get_elem_solution_derivative(), 1.0); for (unsigned int j=0; j != n_u_dofs; j++) { @@ -449,6 +462,10 @@ { FEMContext &c = libmesh_cast_ref(context); + FEBase* p_elem_fe; + + c.get_element_fe( p_var, p_elem_fe ); + // Pin p = 0 at the origin const Point zero(0.,0.); @@ -457,15 +474,15 @@ // The pressure penalty value. \f$ \frac{1}{\epsilon} \f$ const Real penalty = 1.e9; - DenseSubMatrix &Kpp = *c.elem_subjacobians[p_var][p_var]; - DenseSubVector &Fp = *c.elem_subresiduals[p_var]; - const unsigned int n_p_dofs = c.dof_indices_var[p_var].size(); + DenseSubMatrix &Kpp = c.get_elem_jacobian( p_var, p_var ); + DenseSubVector &Fp = c.get_elem_residual( p_var ); + const unsigned int n_p_dofs = c.get_dof_indices( p_var ).size(); Number p = c.point_value(p_var, zero); Number p_value = 0.; unsigned int dim = get_mesh().mesh_dimension(); - FEType fe_type = c.element_fe_var[p_var]->get_fe_type(); + FEType fe_type = p_elem_fe->get_fe_type(); Point p_master = FEInterface::inverse_map(dim, fe_type, c.elem, zero); std::vector point_phi(n_p_dofs); @@ -477,9 +494,9 @@ for (unsigned int i=0; i != n_p_dofs; i++) { Fp(i) += penalty * (p - p_value) * point_phi[i]; - if (request_jacobian && c.elem_solution_derivative) + if (request_jacobian && c.get_elem_solution_derivative()) { - libmesh_assert_equal_to (c.elem_solution_derivative, 1.0); + libmesh_assert_equal_to (c.get_elem_solution_derivative(), 1.0); for (unsigned int j=0; j != n_p_dofs; j++) Kpp(i,j) += penalty * point_phi[i] * point_phi[j]; @@ -503,26 +520,29 @@ { FEMContext &c = libmesh_cast_ref(context); + FEBase* u_elem_fe; + + c.get_element_fe( u_var, u_elem_fe ); + // The subvectors and submatrices we need to fill: const unsigned int dim = this->get_mesh().mesh_dimension(); // Element Jacobian * quadrature weight for interior integration - const std::vector &JxW = c.element_fe_var[u_var]->get_JxW(); + const std::vector &JxW = u_elem_fe->get_JxW(); // The velocity shape functions at interior quadrature points. - const std::vector >& phi = - c.element_fe_var[u_var]->get_phi(); + const std::vector >& phi = u_elem_fe->get_phi(); // The subvectors and submatrices we need to fill: - DenseSubVector &Fu = *c.elem_subresiduals[u_var]; - DenseSubVector &Fv = *c.elem_subresiduals[v_var]; - DenseSubVector &Fw = *c.elem_subresiduals[w_var]; - DenseSubMatrix &Kuu = *c.elem_subjacobians[u_var][u_var]; - DenseSubMatrix &Kvv = *c.elem_subjacobians[v_var][v_var]; - DenseSubMatrix &Kww = *c.elem_subjacobians[w_var][w_var]; + DenseSubVector &Fu = c.get_elem_residual( u_var ); + DenseSubVector &Fv = c.get_elem_residual( v_var ); + DenseSubVector &Fw = c.get_elem_residual( w_var ); + DenseSubMatrix &Kuu = c.get_elem_jacobian( u_var, u_var ); + DenseSubMatrix &Kvv = c.get_elem_jacobian( v_var, v_var ); + DenseSubMatrix &Kww = c.get_elem_jacobian( w_var, w_var ); // The number of local degrees of freedom in velocity - const unsigned int n_u_dofs = c.dof_indices_var[u_var].size(); + const unsigned int n_u_dofs = c.get_dof_indices( u_var ).size(); unsigned int n_qpoints = (c.get_element_qrule())->n_points(); @@ -545,9 +565,9 @@ if (dim == 3) Fw(i) += JxWxRexW * phi[i][qp]; - if (request_jacobian && c.elem_solution_derivative) + if (request_jacobian && c.get_elem_solution_derivative()) { - libmesh_assert_equal_to (c.elem_solution_derivative, 1.0); + libmesh_assert_equal_to (c.get_elem_solution_derivative(), 1.0); Number JxWxRexPhiI = JxWxRe * phi[i][qp]; Number JxWxRexPhiII = JxWxRexPhiI * phi[i][qp]; Index: include/systems/diff_context.h =================================================================== --- include/systems/diff_context.h (revision 6229) +++ include/systems/diff_context.h (working copy) @@ -85,6 +85,81 @@ */ virtual void elem_edge_reinit(Real) {} + const DenseVector& get_elem_solution() const + { return elem_solution; } + + const DenseSubVector& get_elem_solution( unsigned int var ) const + { return *(elem_subsolutions[var]); } + + const DenseVector& get_elem_fixed_solution() const + { return elem_fixed_solution; } + + const DenseSubVector& get_elem_fixed_solution( unsigned int var ) const + { return *(elem_fixed_subsolutions[var]); } + + const DenseVector& get_elem_residual() const + { return elem_residual; } + + DenseVector& get_elem_residual() + { return elem_residual; } + + const DenseSubVector& get_elem_residual( unsigned int var ) const + { return *(elem_subresiduals[var]); } + + DenseSubVector& get_elem_residual( unsigned int var ) + { return *(elem_subresiduals[var]); } + + const DenseMatrix& get_elem_jacobian() const + { return elem_jacobian; } + + DenseMatrix& get_elem_jacobian() + { return elem_jacobian; } + + const DenseSubMatrix& get_elem_jacobian( unsigned int var1, unsigned int var2 ) const + { return *(elem_subjacobians[var1][var2]); } + + DenseSubMatrix& get_elem_jacobian( unsigned int var1, unsigned int var2 ) + { return *(elem_subjacobians[var1][var2]); } + + const std::vector& get_qois() const + { return elem_qoi; } + + std::vector& get_qois() + { return elem_qoi; } + + const std::vector > & get_qoi_derivatives() const + { return elem_qoi_derivative; } + + std::vector > & get_qoi_derivatives() + { return elem_qoi_derivative; } + + const DenseSubVector& get_qoi_derivatives( unsigned int qoi, unsigned int var ) const + { return *(elem_qoi_subderivatives[qoi][var]); } + + DenseSubVector& get_qoi_derivatives( unsigned int qoi, unsigned int var ) + { return *(elem_qoi_subderivatives[qoi][var]); } + + const std::vector& get_dof_indices() const + { return dof_indices; } + + const std::vector& get_dof_indices( unsigned int var ) const + { return dof_indices_var[var]; } + + Real get_system_time() const + { return system_time; } + + Real get_time() const + { return time; } + + void set_time( Real time_in ) + { time = time_in; } + + Real get_elem_solution_derivative() const + { return elem_solution_derivative; } + + Real get_fixed_solution_derivative() const + { return fixed_solution_derivative; } + /** * For time-dependent problems, this is the time t for which the current * nonlinear_solution is defined.