Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
l2.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_integrators_l2_h
17 #define dealii_integrators_l2_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/quadrature.h>
24 
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/fe/mapping.h>
27 
28 #include <deal.II/lac/full_matrix.h>
29 
30 #include <deal.II/meshworker/dof_info.h>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 namespace LocalIntegrators
35 {
43  namespace L2
44  {
61  template <int dim>
62  void
64  const FEValuesBase<dim> &fe,
65  const double factor = 1.)
66  {
67  const unsigned int n_dofs = fe.dofs_per_cell;
68  const unsigned int n_components = fe.get_fe().n_components();
69 
70  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
71  {
72  const double dx = fe.JxW(k) * factor;
73  for (unsigned int i = 0; i < n_dofs; ++i)
74  {
75  double Mii = 0.0;
76  for (unsigned int d = 0; d < n_components; ++d)
77  Mii += dx * fe.shape_value_component(i, k, d) *
78  fe.shape_value_component(i, k, d);
79 
80  M(i, i) += Mii;
81 
82  for (unsigned int j = i + 1; j < n_dofs; ++j)
83  {
84  double Mij = 0.0;
85  for (unsigned int d = 0; d < n_components; ++d)
86  Mij += dx * fe.shape_value_component(j, k, d) *
87  fe.shape_value_component(i, k, d);
88 
89  M(i, j) += Mij;
90  M(j, i) += Mij;
91  }
92  }
93  }
94  }
95 
115  template <int dim>
116  void
118  const FEValuesBase<dim> & fe,
119  const std::vector<double> &weights)
120  {
121  const unsigned int n_dofs = fe.dofs_per_cell;
122  const unsigned int n_components = fe.get_fe().n_components();
123  AssertDimension(M.m(), n_dofs);
124  AssertDimension(M.n(), n_dofs);
125  AssertDimension(weights.size(), fe.n_quadrature_points);
126 
127  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
128  {
129  const double dx = fe.JxW(k) * weights[k];
130  for (unsigned int i = 0; i < n_dofs; ++i)
131  {
132  double Mii = 0.0;
133  for (unsigned int d = 0; d < n_components; ++d)
134  Mii += dx * fe.shape_value_component(i, k, d) *
135  fe.shape_value_component(i, k, d);
136 
137  M(i, i) += Mii;
138 
139  for (unsigned int j = i + 1; j < n_dofs; ++j)
140  {
141  double Mij = 0.0;
142  for (unsigned int d = 0; d < n_components; ++d)
143  Mij += dx * fe.shape_value_component(j, k, d) *
144  fe.shape_value_component(i, k, d);
145 
146  M(i, j) += Mij;
147  M(j, i) += Mij;
148  }
149  }
150  }
151  }
152 
169  template <int dim, typename number>
170  void
171  L2(Vector<number> & result,
172  const FEValuesBase<dim> & fe,
173  const std::vector<double> &input,
174  const double factor = 1.)
175  {
176  const unsigned int n_dofs = fe.dofs_per_cell;
177  AssertDimension(result.size(), n_dofs);
178  AssertDimension(fe.get_fe().n_components(), 1);
179  AssertDimension(input.size(), fe.n_quadrature_points);
180 
181  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
182  for (unsigned int i = 0; i < n_dofs; ++i)
183  result(i) += fe.JxW(k) * factor * input[k] * fe.shape_value(i, k);
184  }
185 
202  template <int dim, typename number>
203  void
204  L2(Vector<number> & result,
205  const FEValuesBase<dim> & fe,
206  const ArrayView<const std::vector<double>> &input,
207  const double factor = 1.)
208  {
209  const unsigned int n_dofs = fe.dofs_per_cell;
210  const unsigned int n_components = input.size();
211 
212  AssertDimension(result.size(), n_dofs);
213  AssertDimension(input.size(), fe.get_fe().n_components());
214 
215  for (unsigned int k = 0; k < fe.n_quadrature_points; ++k)
216  for (unsigned int i = 0; i < n_dofs; ++i)
217  for (unsigned int d = 0; d < n_components; ++d)
218  result(i) += fe.JxW(k) * factor *
219  fe.shape_value_component(i, k, d) * input[d][k];
220  }
221 
253  template <int dim>
254  void
256  FullMatrix<double> & M12,
257  FullMatrix<double> & M21,
258  FullMatrix<double> & M22,
259  const FEValuesBase<dim> &fe1,
260  const FEValuesBase<dim> &fe2,
261  const double factor1 = 1.,
262  const double factor2 = 1.)
263  {
264  const unsigned int n1_dofs = fe1.dofs_per_cell;
265  const unsigned int n2_dofs = fe2.dofs_per_cell;
266  const unsigned int n_components = fe1.get_fe().n_components();
267 
268  Assert(n1_dofs == n2_dofs, ExcNotImplemented());
269  (void)n2_dofs;
270  AssertDimension(n_components, fe2.get_fe().n_components());
271  AssertDimension(M11.m(), n1_dofs);
272  AssertDimension(M12.m(), n1_dofs);
273  AssertDimension(M21.m(), n2_dofs);
274  AssertDimension(M22.m(), n2_dofs);
275  AssertDimension(M11.n(), n1_dofs);
276  AssertDimension(M12.n(), n2_dofs);
277  AssertDimension(M21.n(), n1_dofs);
278  AssertDimension(M22.n(), n2_dofs);
279 
280  for (unsigned int k = 0; k < fe1.n_quadrature_points; ++k)
281  {
282  const double dx = fe1.JxW(k);
283 
284  for (unsigned int i = 0; i < n1_dofs; ++i)
285  for (unsigned int j = 0; j < n1_dofs; ++j)
286  for (unsigned int d = 0; d < n_components; ++d)
287  {
288  const double u1 =
289  factor1 * fe1.shape_value_component(j, k, d);
290  const double u2 =
291  -factor2 * fe2.shape_value_component(j, k, d);
292  const double v1 =
293  factor1 * fe1.shape_value_component(i, k, d);
294  const double v2 =
295  -factor2 * fe2.shape_value_component(i, k, d);
296 
297  M11(i, j) += dx * u1 * v1;
298  M12(i, j) += dx * u2 * v1;
299  M21(i, j) += dx * u1 * v2;
300  M22(i, j) += dx * u2 * v2;
301  }
302  }
303  }
304  } // namespace L2
305 } // namespace LocalIntegrators
306 
307 DEAL_II_NAMESPACE_CLOSE
308 
309 #endif
size_type m() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
const unsigned int dofs_per_cell
Definition: fe_values.h:2106
const FiniteElement< dim, spacedim > & get_fe() const
void jump_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const double factor1=1., const double factor2=1.)
Definition: l2.h:255
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
Library of integrals over cells and faces.
Definition: advection.h:34
void weighted_mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const std::vector< double > &weights)
Definition: l2.h:117
size_type n() const
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
#define Assert(cond, exc)
Definition: exceptions.h:1407
const unsigned int n_quadrature_points
Definition: fe_values.h:2099
double JxW(const unsigned int quadrature_point) const
size_type size() const
static ::ExceptionBase & ExcNotImplemented()
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:171