MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mortarintegrator.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "mortarintegrator.hpp"
13 
14 namespace mfem
15 {
16 
18  const FiniteElement &trial, const IntegrationRule &trial_ir,
19  ElementTransformation &trial_Trans, const FiniteElement &test,
20  const IntegrationRule &test_ir, ElementTransformation &test_Trans,
21  DenseMatrix &elmat)
22 {
23  int tr_nd = trial.GetDof();
24  int te_nd = test.GetDof();
25  double w;
26 
27  Vector shape, te_shape;
28 
29  elmat.SetSize(te_nd, tr_nd);
30  shape.SetSize(tr_nd);
31  te_shape.SetSize(te_nd);
32 
33  elmat = 0.0;
34  for (int i = 0; i < test_ir.GetNPoints(); i++)
35  {
36  const IntegrationPoint &trial_ip = trial_ir.IntPoint(i);
37  const IntegrationPoint &test_ip = test_ir.IntPoint(i);
38  test_Trans.SetIntPoint(&test_ip);
39 
40  trial.CalcShape(trial_ip, shape);
41  test.CalcShape(test_ip, te_shape);
42 
43  w = test_Trans.Weight() * test_ip.weight;
44 
45  te_shape *= w;
46  AddMultVWt(te_shape, shape, elmat);
47  }
48 }
49 
51  const FiniteElement &trial, const IntegrationRule &trial_ir,
52  ElementTransformation &trial_Trans, const FiniteElement &test,
53  const IntegrationRule &test_ir, ElementTransformation &test_Trans,
54  DenseMatrix &elmat)
55 {
56  if (test.GetRangeType() == FiniteElement::SCALAR && VQ)
57  {
58  // assume test is scalar FE and trial is vector FE
59  int dim = test.GetDim();
60  int trial_dof = trial.GetDof();
61  int test_dof = test.GetDof();
62  double w;
63 
64  if (MQ)
65  mfem_error("VectorFEMassIntegrator::AssembleElementMatrix2(...)\n"
66  " is not implemented for tensor materials");
67 
68 #ifdef MFEM_THREAD_SAFE
69  DenseMatrix trial_vshape(trial_dof, dim);
70  Vector shape(test_dof);
71  Vector D(dim);
72 #else
73  trial_vshape.SetSize(trial_dof, dim);
74  shape.SetSize(test_dof);
75  D.SetSize(dim);
76 #endif
77 
78  elmat.SetSize(test_dof, trial_dof);
79 
80  elmat = 0.0;
81  for (int i = 0; i < test_ir.GetNPoints(); i++)
82  {
83  const IntegrationPoint &trial_ip = trial_ir.IntPoint(i);
84  const IntegrationPoint &test_ip = test_ir.IntPoint(i);
85 
86  trial_Trans.SetIntPoint(&trial_ip);
87  test_Trans.SetIntPoint(&test_ip);
88 
89  trial.CalcVShape(trial_Trans, trial_vshape);
90  test.CalcShape(test_ip, shape);
91 
92  w = test_ip.weight * test_Trans.Weight();
93  VQ->Eval(D, test_Trans, test_ip);
94  D *= w;
95 
96  for (int d = 0; d < dim; d++)
97  {
98  for (int j = 0; j < test_dof; j++)
99  {
100  for (int k = 0; k < trial_dof; k++)
101  {
102  elmat(j, k) += D[d] * shape(j) * trial_vshape(k, d);
103  }
104  }
105  }
106  }
107  }
108  else if (test.GetRangeType() == FiniteElement::SCALAR)
109  {
110  // assume test is scalar FE and trial is vector FE
111  int dim = test.GetDim();
112  int trial_dof = trial.GetDof();
113  int test_dof = test.GetDof();
114  double w;
115 
116  if (VQ || MQ)
117  mfem_error("VectorFEMassIntegrator::AssembleElementMatrix2(...)\n"
118  " is not implemented for vector/tensor permeability");
119 
120 #ifdef MFEM_THREAD_SAFE
121  DenseMatrix trial_vshape(trial_dof, dim);
122  Vector shape(test_dof);
123 #else
124  trial_vshape.SetSize(trial_dof, dim);
125  shape.SetSize(test_dof);
126 #endif
127 
128  elmat.SetSize(dim * test_dof, trial_dof);
129 
130  elmat = 0.0;
131  for (int i = 0; i < test_ir.GetNPoints(); i++)
132  {
133  const IntegrationPoint &trial_ip = trial_ir.IntPoint(i);
134  const IntegrationPoint &test_ip = test_ir.IntPoint(i);
135 
136  trial_Trans.SetIntPoint(&trial_ip);
137  test_Trans.SetIntPoint(&test_ip);
138 
139  trial.CalcVShape(trial_Trans, trial_vshape);
140  test.CalcShape(test_ip, shape);
141 
142  w = test_ip.weight * test_Trans.Weight();
143 
144  if (Q)
145  {
146  w *= Q->Eval(test_Trans, test_ip);
147  }
148 
149  for (int d = 0; d < dim; d++)
150  {
151  for (int j = 0; j < test_dof; j++)
152  {
153  for (int k = 0; k < trial_dof; k++)
154  {
155  elmat(d * test_dof + j, k) += w * shape(j) * trial_vshape(k, d);
156  }
157  }
158  }
159  }
160  }
161  else
162  {
163  // assume both test and trial are vector FE
164  int dim = test.GetDim();
165  int trial_dof = trial.GetDof();
166  int test_dof = test.GetDof();
167  double w;
168 
169  if (VQ || MQ)
170  mfem_error("VectorFEMassIntegrator::AssembleElementMatrix2(...)\n"
171  " is not implemented for vector/tensor permeability");
172 
173 #ifdef MFEM_THREAD_SAFE
174  DenseMatrix trial_vshape(trial_dof, dim);
175  DenseMatrix test_vshape(test_dof, dim);
176 #else
177  trial_vshape.SetSize(trial_dof, dim);
178  test_vshape.SetSize(test_dof, dim);
179 #endif
180 
181  elmat.SetSize(test_dof, trial_dof);
182 
183  elmat = 0.0;
184  for (int i = 0; i < test_ir.GetNPoints(); i++)
185  {
186  const IntegrationPoint &trial_ip = trial_ir.IntPoint(i);
187  const IntegrationPoint &test_ip = test_ir.IntPoint(i);
188 
189  trial_Trans.SetIntPoint(&trial_ip);
190  test_Trans.SetIntPoint(&test_ip);
191 
192  trial.CalcVShape(trial_Trans, trial_vshape);
193  test.CalcVShape(test_Trans, test_vshape);
194 
195  w = test_ip.weight * test_Trans.Weight();
196  if (Q)
197  {
198  w *= Q->Eval(test_Trans, test_ip);
199  }
200 
201  for (int d = 0; d < dim; d++)
202  {
203  for (int j = 0; j < test_dof; j++)
204  {
205  for (int k = 0; k < trial_dof; k++)
206  {
207  elmat(j, k) += w * test_vshape(j, d) * trial_vshape(k, d);
208  }
209  }
210  }
211  }
212  }
213 }
214 
215 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3127
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_base.cpp:39
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void AssembleElementMatrix(const FiniteElement &trial, const IntegrationRule &trial_ir, ElementTransformation &trial_Trans, const FiniteElement &test, const IntegrationRule &test_ir, ElementTransformation &test_Trans, DenseMatrix &elemmat) override
Implements the assembly routine.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Vector data type.
Definition: vector.hpp:60
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:340
void AssembleElementMatrix(const FiniteElement &trial, const IntegrationRule &trial_ir, ElementTransformation &trial_Trans, const FiniteElement &test, const IntegrationRule &test_ir, ElementTransformation &test_Trans, DenseMatrix &elemmat) override
Implements the assembly routine.