MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mortarintegrator.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
14namespace 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
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
Abstract class for all finite elements.
Definition fe_base.hpp:239
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:40
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:316
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition fe_base.hpp:346
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...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:329
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
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.
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 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.
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
int dim
Definition ex24.cpp:53
void mfem_error(const char *msg)
Definition error.cpp:154
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.