Reference documentation for deal.II version 8.1.0
laplace.h
1 // ---------------------------------------------------------------------
2 // @f$Id: laplace.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2010 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__integrators_laplace_h
18 #define __deal2__integrators_laplace_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/fe/mapping.h>
26 #include <deal.II/fe/fe_values.h>
27 #include <deal.II/meshworker/dof_info.h>
28 
30 
31 namespace LocalIntegrators
32 {
40  namespace Laplace
41  {
55  template<int dim>
56  void cell_matrix (
58  const FEValuesBase<dim> &fe,
59  const double factor = 1.)
60  {
61  const unsigned int n_dofs = fe.dofs_per_cell;
62  const unsigned int n_components = fe.get_fe().n_components();
63 
64  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
65  {
66  const double dx = fe.JxW(k) * factor;
67  for (unsigned int i=0; i<n_dofs; ++i)
68  {
69  for (unsigned int j=0; j<n_dofs; ++j)
70  for (unsigned int d=0; d<n_components; ++d)
71  M(i,j) += dx *
72  (fe.shape_grad_component(j,k,d) * fe.shape_grad_component(i,k,d));
73  }
74  }
75  }
76 
84  template <int dim>
85  inline void
87  Vector<double> &result,
88  const FEValuesBase<dim> &fe,
89  const std::vector<Tensor<1,dim> > &input,
90  double factor = 1.)
91  {
92  const unsigned int nq = fe.n_quadrature_points;
93  const unsigned int n_dofs = fe.dofs_per_cell;
94  Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
95  Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
96 
97  for (unsigned int k=0; k<nq; ++k)
98  {
99  const double dx = factor * fe.JxW(k);
100  for (unsigned int i=0; i<n_dofs; ++i)
101  result(i) += dx * (input[k] * fe.shape_grad(i,k));
102  }
103  }
104 
105 
113  template <int dim>
114  inline void
116  Vector<double> &result,
117  const FEValuesBase<dim> &fe,
118  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &input,
119  double factor = 1.)
120  {
121  const unsigned int nq = fe.n_quadrature_points;
122  const unsigned int n_dofs = fe.dofs_per_cell;
123  const unsigned int n_comp = fe.get_fe().n_components();
124 
126  Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
127 
128  for (unsigned int k=0; k<nq; ++k)
129  {
130  const double dx = factor * fe.JxW(k);
131  for (unsigned int i=0; i<n_dofs; ++i)
132  for (unsigned int d=0; d<n_comp; ++d)
133  {
134 
135  result(i) += dx * (input[d][k] * fe.shape_grad_component(i,k,d));
136  }
137  }
138  }
139 
140 
153  template <int dim>
156  const FEValuesBase<dim> &fe,
157  double penalty,
158  double factor = 1.)
159  {
160  const unsigned int n_dofs = fe.dofs_per_cell;
161  const unsigned int n_comp = fe.get_fe().n_components();
162 
163  Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs));
164  Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs));
165 
166  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
167  {
168  const double dx = fe.JxW(k) * factor;
169  const Point<dim> &n = fe.normal_vector(k);
170  for (unsigned int i=0; i<n_dofs; ++i)
171  for (unsigned int j=0; j<n_dofs; ++j)
172  for (unsigned int d=0; d<n_comp; ++d)
173  M(i,j) += dx *
174  (2. * fe.shape_value_component(i,k,d) * penalty * fe.shape_value_component(j,k,d)
175  - (n * fe.shape_grad_component(i,k,d)) * fe.shape_value_component(j,k,d)
176  - (n * fe.shape_grad_component(j,k,d)) * fe.shape_value_component(i,k,d));
177  }
178  }
179 
196  template <int dim>
198  Vector<double> &result,
199  const FEValuesBase<dim> &fe,
200  const std::vector<double> &input,
201  const std::vector<Tensor<1,dim> > &Dinput,
202  const std::vector<double> &data,
203  double penalty,
204  double factor = 1.)
205  {
206  const unsigned int n_dofs = fe.dofs_per_cell;
207  AssertDimension(input.size(), fe.n_quadrature_points);
208  AssertDimension(Dinput.size(), fe.n_quadrature_points);
209  AssertDimension(data.size(), fe.n_quadrature_points);
210 
211  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
212  {
213  const double dx = factor * fe.JxW(k);
214  const Point<dim> &n = fe.normal_vector(k);
215  for (unsigned int i=0; i<n_dofs; ++i)
216  {
217  const double dnv = fe.shape_grad(i,k) * n;
218  const double dnu = Dinput[k] * n;
219  const double v= fe.shape_value(i,k);
220  const double u= input[k];
221  const double g= data[k];
222 
223  result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
224  }
225  }
226  }
227 
247  template <int dim>
249  Vector<double> &result,
250  const FEValuesBase<dim> &fe,
251  const VectorSlice<const std::vector<std::vector<double> > > &input,
252  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput,
253  const VectorSlice<const std::vector<std::vector<double> > > &data,
254  double penalty,
255  double factor = 1.)
256  {
257  const unsigned int n_dofs = fe.dofs_per_cell;
258  const unsigned int n_comp = fe.get_fe().n_components();
262 
263  for (unsigned int k=0; k<fe.n_quadrature_points; ++k)
264  {
265  const double dx = factor * fe.JxW(k);
266  const Point<dim> &n = fe.normal_vector(k);
267  for (unsigned int i=0; i<n_dofs; ++i)
268  for (unsigned int d=0; d<n_comp; ++d)
269  {
270  const double dnv = fe.shape_grad_component(i,k,d) * n;
271  const double dnu = Dinput[d][k] * n;
272  const double v= fe.shape_value_component(i,k,d);
273  const double u= input[d][k];
274  const double g= data[d][k];
275 
276  result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
277  }
278  }
279  }
280 
300  template <int dim>
301  void ip_matrix (
302  FullMatrix<double> &M11,
303  FullMatrix<double> &M12,
304  FullMatrix<double> &M21,
305  FullMatrix<double> &M22,
306  const FEValuesBase<dim> &fe1,
307  const FEValuesBase<dim> &fe2,
308  double penalty,
309  double factor1 = 1.,
310  double factor2 = -1.)
311  {
312  const unsigned int n_dofs = fe1.dofs_per_cell;
313  AssertDimension(M11.n(), n_dofs);
314  AssertDimension(M11.m(), n_dofs);
315  AssertDimension(M12.n(), n_dofs);
316  AssertDimension(M12.m(), n_dofs);
317  AssertDimension(M21.n(), n_dofs);
318  AssertDimension(M21.m(), n_dofs);
319  AssertDimension(M22.n(), n_dofs);
320  AssertDimension(M22.m(), n_dofs);
321 
322  const double nui = factor1;
323  const double nue = (factor2 < 0) ? factor1 : factor2;
324  const double nu = .5*(nui+nue);
325 
326  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
327  {
328  const double dx = fe1.JxW(k);
329  const Point<dim> &n = fe1.normal_vector(k);
330  for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
331  {
332  for (unsigned int i=0; i<n_dofs; ++i)
333  {
334  for (unsigned int j=0; j<n_dofs; ++j)
335  {
336  const double vi = fe1.shape_value_component(i,k,d);
337  const double dnvi = n * fe1.shape_grad_component(i,k,d);
338  const double ve = fe2.shape_value_component(i,k,d);
339  const double dnve = n * fe2.shape_grad_component(i,k,d);
340  const double ui = fe1.shape_value_component(j,k,d);
341  const double dnui = n * fe1.shape_grad_component(j,k,d);
342  const double ue = fe2.shape_value_component(j,k,d);
343  const double dnue = n * fe2.shape_grad_component(j,k,d);
344  M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
345  M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
346  M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
347  M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
348  }
349  }
350  }
351  }
352  }
353 
368  template <int dim>
370  FullMatrix<double> &M11,
371  FullMatrix<double> &M12,
372  FullMatrix<double> &M21,
373  FullMatrix<double> &M22,
374  const FEValuesBase<dim> &fe1,
375  const FEValuesBase<dim> &fe2,
376  double penalty,
377  double factor1 = 1.,
378  double factor2 = -1.)
379  {
380  const unsigned int n_dofs = fe1.dofs_per_cell;
381  AssertDimension(fe1.get_fe().n_components(), dim);
382  AssertDimension(fe2.get_fe().n_components(), dim);
383  AssertDimension(M11.n(), n_dofs);
384  AssertDimension(M11.m(), n_dofs);
385  AssertDimension(M12.n(), n_dofs);
386  AssertDimension(M12.m(), n_dofs);
387  AssertDimension(M21.n(), n_dofs);
388  AssertDimension(M21.m(), n_dofs);
389  AssertDimension(M22.n(), n_dofs);
390  AssertDimension(M22.m(), n_dofs);
391 
392  const double nui = factor1;
393  const double nue = (factor2 < 0) ? factor1 : factor2;
394  const double nu = .5*(nui+nue);
395 
396  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
397  {
398  const double dx = fe1.JxW(k);
399  const Point<dim> &n = fe1.normal_vector(k);
400  for (unsigned int i=0; i<n_dofs; ++i)
401  {
402  for (unsigned int j=0; j<n_dofs; ++j)
403  {
404  double u1dotn = 0.;
405  double v1dotn = 0.;
406  double u2dotn = 0.;
407  double v2dotn = 0.;
408 
409  double ngradu1n = 0.;
410  double ngradv1n = 0.;
411  double ngradu2n = 0.;
412  double ngradv2n = 0.;
413 
414  for (unsigned int d=0; d<dim; ++d)
415  {
416  u1dotn += n(d)*fe1.shape_value_component(j,k,d);
417  v1dotn += n(d)*fe1.shape_value_component(i,k,d);
418  u2dotn += n(d)*fe2.shape_value_component(j,k,d);
419  v2dotn += n(d)*fe2.shape_value_component(i,k,d);
420 
421  ngradu1n += n*fe1.shape_grad_component(j,k,d)*n(d);
422  ngradv1n += n*fe1.shape_grad_component(i,k,d)*n(d);
423  ngradu2n += n*fe2.shape_grad_component(j,k,d)*n(d);
424  ngradv2n += n*fe2.shape_grad_component(i,k,d)*n(d);
425  }
426 
427  for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
428  {
429  const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n(d);
430  const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n(d);
431 
432  const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n(d);
433  const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n(d);
434 
435  const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n(d);
436  const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n(d);
437 
438  const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n(d);
439  const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n(d);
440 
441  M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
442  M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
443  M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
444  M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
445  }
446  }
447  }
448  }
449  }
450 
461  template<int dim>
462  void
464  Vector<double> &result1,
465  Vector<double> &result2,
466  const FEValuesBase<dim> &fe1,
467  const FEValuesBase<dim> &fe2,
468  const std::vector<double> &input1,
469  const std::vector<Tensor<1,dim> > &Dinput1,
470  const std::vector<double> &input2,
471  const std::vector<Tensor<1,dim> > &Dinput2,
472  double pen,
473  double int_factor = 1.,
474  double ext_factor = -1.)
475  {
476  Assert(fe1.get_fe().n_components() == 1,
477  ExcDimensionMismatch(fe1.get_fe().n_components(), 1));
478  Assert(fe2.get_fe().n_components() == 1,
479  ExcDimensionMismatch(fe2.get_fe().n_components(), 1));
480 
481  const double nui = int_factor;
482  const double nue = (ext_factor < 0) ? int_factor : ext_factor;
483  const double penalty = .5 * pen * (nui + nue);
484 
485  const unsigned int n_dofs = fe1.dofs_per_cell;
486 
487  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
488  {
489  const double dx = fe1.JxW(k);
490  const Point<dim> &n = fe1.normal_vector(k);
491 
492  for (unsigned int i=0; i<n_dofs; ++i)
493  {
494  const double vi = fe1.shape_value(i,k);
495  const Tensor<1,dim> &Dvi = fe1.shape_grad(i,k);
496  const double dnvi = Dvi * n;
497  const double ve = fe2.shape_value(i,k);
498  const Tensor<1,dim> &Dve = fe2.shape_grad(i,k);
499  const double dnve = Dve * n;
500 
501  const double ui = input1[k];
502  const Tensor<1,dim> &Dui = Dinput1[k];
503  const double dnui = Dui * n;
504  const double ue = input2[k];
505  const Tensor<1,dim> &Due = Dinput2[k];
506  const double dnue = Due * n;
507 
508  result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
509  result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
510  result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
511  result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
512  }
513  }
514  }
515 
516 
528  template<int dim>
529  void
531  Vector<double> &result1,
532  Vector<double> &result2,
533  const FEValuesBase<dim> &fe1,
534  const FEValuesBase<dim> &fe2,
535  const VectorSlice<const std::vector<std::vector<double> > > &input1,
536  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput1,
537  const VectorSlice<const std::vector<std::vector<double> > > &input2,
538  const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > > &Dinput2,
539  double pen,
540  double int_factor = 1.,
541  double ext_factor = -1.)
542  {
543  const unsigned int n_comp = fe1.get_fe().n_components();
544  const unsigned int n1 = fe1.dofs_per_cell;
545 
547  AssertVectorVectorDimension(Dinput1, n_comp, fe1.n_quadrature_points);
549  AssertVectorVectorDimension(Dinput2, n_comp, fe2.n_quadrature_points);
550 
551  const double nui = int_factor;
552  const double nue = (ext_factor < 0) ? int_factor : ext_factor;
553  const double penalty = .5 * pen * (nui + nue);
554 
555 
556  for (unsigned int k=0; k<fe1.n_quadrature_points; ++k)
557  {
558  const double dx = fe1.JxW(k);
559  const Point<dim> &n = fe1.normal_vector(k);
560 
561  for (unsigned int i=0; i<n1; ++i)
562  for (unsigned int d=0; d<n_comp; ++d)
563  {
564  const double vi = fe1.shape_value_component(i,k,d);
565  const Tensor<1,dim> &Dvi = fe1.shape_grad_component(i,k,d);
566  const double dnvi = Dvi * n;
567  const double ve = fe2.shape_value_component(i,k,d);
568  const Tensor<1,dim> &Dve = fe2.shape_grad_component(i,k,d);
569  const double dnve = Dve * n;
570 
571  const double ui = input1[d][k];
572  const Tensor<1,dim> &Dui = Dinput1[d][k];
573  const double dnui = Dui * n;
574  const double ue = input2[d][k];
575  const Tensor<1,dim> &Due = Dinput2[d][k];
576  const double dnue = Due * n;
577 
578  result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
579  result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
580  result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
581  result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
582  }
583  }
584  }
585 
586 
587 
602  template <int dim, int spacedim, typename number>
606  unsigned int deg1,
607  unsigned int deg2)
608  {
609  const unsigned int normal1 = GeometryInfo<dim>::unit_normal_direction[dinfo1.face_number];
610  const unsigned int normal2 = GeometryInfo<dim>::unit_normal_direction[dinfo2.face_number];
611  const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1+1);
612  const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2+1);
613 
614  double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
615  double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
616  if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children())
617  {
618  Assert (dinfo1.face == dinfo2.face, ExcInternalError());
619  Assert (dinfo1.face->has_children(), ExcInternalError());
620  penalty1 *= 2;
621  }
622  const double penalty = 0.5*(penalty1 + penalty2);
623  return penalty;
624  }
625  }
626 }
627 
628 
629 DEAL_II_NAMESPACE_CLOSE
630 
631 #endif
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
const Tensor< 1, spacedim > & shape_grad(const unsigned int function_no, const unsigned int quadrature_point) const
Triangulation< dim, spacedim >::cell_iterator cell
The current cell.
Definition: dof_info.h:73
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
Tensor< 1, spacedim > shape_grad_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
void ip_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, double penalty, double factor1=1., double factor2=-1.)
Definition: laplace.h:301
const unsigned int dofs_per_cell
Definition: fe_values.h:1418
const Point< spacedim > & normal_vector(const unsigned int i) const
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: laplace.h:56
#define AssertVectorVectorDimension(vec, dim1, dim2)
Definition: exceptions.h:870
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, double factor=1.)
Definition: laplace.h:86
Library of integrals over cells and faces.
Definition: elasticity.h:31
void nitsche_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const std::vector< Tensor< 1, dim > > &Dinput, const std::vector< double > &data, double penalty, double factor=1.)
Definition: laplace.h:197
void ip_residual(Vector< double > &result1, Vector< double > &result2, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, const std::vector< double > &input1, const std::vector< Tensor< 1, dim > > &Dinput1, const std::vector< double > &input2, const std::vector< Tensor< 1, dim > > &Dinput2, double pen, double int_factor=1., double ext_factor=-1.)
Definition: laplace.h:463
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t size() const
size_type m() const
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
void nitsche_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, double penalty, double factor=1.)
Definition: laplace.h:154
double compute_penalty(const MeshWorker::DoFInfo< dim, spacedim, number > &dinfo1, const MeshWorker::DoFInfo< dim, spacedim, number > &dinfo2, unsigned int deg1, unsigned int deg2)
Definition: laplace.h:603
unsigned int face_number
Definition: dof_info.h:88
const FiniteElement< dim, spacedim > & get_fe() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcInternalError()
Triangulation< dim, spacedim >::face_iterator face
The current face.
Definition: dof_info.h:76
double JxW(const unsigned int quadrature_point) const
void ip_tangential_matrix(FullMatrix< double > &M11, FullMatrix< double > &M12, FullMatrix< double > &M21, FullMatrix< double > &M22, const FEValuesBase< dim > &fe1, const FEValuesBase< dim > &fe2, double penalty, double factor1=1., double factor2=-1.)
Definition: laplace.h:369