libMesh
L2system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "L2system.h"
19 
20 #include "libmesh/fe_base.h"
21 #include "libmesh/fe_interface.h"
22 #include "libmesh/fem_context.h"
23 #include "libmesh/getpot.h"
24 #include "libmesh/mesh.h"
25 #include "libmesh/quadrature.h"
26 #include "libmesh/string_to_enum.h"
27 
28 using namespace libMesh;
29 
31 {
32  for (std::map<FEMContext *, FEMContext *>::const_iterator
33  it = input_contexts.begin();
34  it != input_contexts.end(); ++it)
35  {
36  FEMContext * c = it->second;
37  delete c;
38  }
39 }
40 
42 {
43  this->add_variable ("u", static_cast<Order>(_fe_order),
44  Utility::string_to_enum<FEFamily>(_fe_family));
45 
46  // Do the parent's initialization after variables are defined
48 }
49 
50 
51 
53 {
54  FEMContext & c = cast_ref<FEMContext &>(context);
55 
56  // Now make sure we have requested all the data
57  // we need to build the L2 system.
58  c.get_element_fe(0)->get_JxW();
59  c.get_element_fe(0)->get_phi();
60  c.get_element_fe(0)->get_xyz();
61 
62  // Build a corresponding context for the input system if we haven't
63  // already
64  FEMContext *& input_context = input_contexts[&c];
65  if (!input_context)
66  {
67  libmesh_assert(input_system);
68  input_context = new FEMContext(*input_system);
69 
70  libmesh_assert(goal_func.get());
71  goal_func->init_context(*input_context);
72  }
73 
74  FEMSystem::init_context(context);
75 }
76 
77 
78 bool L2System::element_time_derivative (bool request_jacobian,
79  DiffContext & context)
80 {
81  FEMContext & c = cast_ref<FEMContext &>(context);
82 
83  // First we get some references to cell-specific data that
84  // will be used to assemble the linear system.
85 
86  // Element Jacobian * quadrature weights for interior integration
87  const std::vector<Real> & JxW = c.get_element_fe(0)->get_JxW();
88 
89  const std::vector<std::vector<Real>> & phi = c.get_element_fe(0)->get_phi();
90 
91  const std::vector<Point> & xyz = c.get_element_fe(0)->get_xyz();
92 
93  // The number of local degrees of freedom in each variable
94  const unsigned int n_u_dofs = c.get_dof_indices(0).size();
95 
96  // The subvectors and submatrices we need to fill:
99 
100  unsigned int n_qpoints = c.get_element_qrule().n_points();
101 
102  libmesh_assert (input_contexts.find(&c) != input_contexts.end());
103 
104  FEMContext & input_c = *input_contexts[&c];
105  input_c.pre_fe_reinit(*input_system, &c.get_elem());
106  input_c.elem_fe_reinit();
107 
108  for (unsigned int qp=0; qp != n_qpoints; qp++)
109  {
110  Number u = c.interior_value(0, qp);
111 
112  Number ufunc = (*goal_func)(input_c, xyz[qp]);
113 
114  for (unsigned int i=0; i != n_u_dofs; i++)
115  F(i) += JxW[qp] * ((u - ufunc) * phi[i][qp]);
116  if (request_jacobian)
117  {
118  const Number JxWxD = JxW[qp] *
120 
121  for (unsigned int i=0; i != n_u_dofs; i++)
122  for (unsigned int j=0; j != n_u_dofs; ++j)
123  K(i,j) += JxWxD * (phi[i][qp] * phi[j][qp]);
124  }
125  } // end of the quadrature point qp-loop
126 
127  return request_jacobian;
128 }
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: L2system.C:41
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1582
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
The libMesh namespace provides an interface to certain functionality in the library.
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
virtual bool element_time_derivative(bool request_jacobian, libMesh::DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: L2system.C:78
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
libmesh_assert(j)
Defines a dense subvector for use in finite element computations.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
Defines a dense submatrix for use in Finite Element-type computations.
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Real get_elem_solution_derivative() const
The derivative of the current elem_solution w.r.t.
Definition: diff_context.h:419
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
~L2System()
Definition: L2system.C:30
unsigned int n_points() const
Definition: quadrature.h:113
virtual void init_context(libMesh::DiffContext &context)
Definition: L2system.C:52
virtual void init_context(DiffContext &) libmesh_override
Definition: fem_system.C:1345