MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_nurbs.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// H1 Finite Element classes utilizing the Bernstein basis
13
14#include "fe_nurbs.hpp"
15#include "../../mesh/nurbs.hpp"
16
17namespace mfem
18{
19
20using namespace std;
21
22void NURBS1DFiniteElement::SetOrder() const
23{
24 order = kv[0]->GetOrder();
25 dof = order + 1;
26
27 weights.SetSize(dof);
28 shape_x.SetSize(dof);
29}
30
31void NURBS1DFiniteElement::CalcShape(const IntegrationPoint &ip,
32 Vector &shape) const
33{
34 kv[0]->CalcShape(shape, ijk[0], ip.x);
35
36 real_t sum = 0.0;
37 for (int i = 0; i <= order; i++)
38 {
39 sum += (shape(i) *= weights(i));
40 }
41
42 shape /= sum;
43}
44
45void NURBS1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
46 DenseMatrix &dshape) const
47{
48 Vector grad(dshape.Data(), dof);
49
50 kv[0]->CalcShape (shape_x, ijk[0], ip.x);
51 kv[0]->CalcDShape(grad, ijk[0], ip.x);
52
53 real_t sum = 0.0, dsum = 0.0;
54 for (int i = 0; i <= order; i++)
55 {
56 sum += (shape_x(i) *= weights(i));
57 dsum += ( grad(i) *= weights(i));
58 }
59
60 sum = 1.0/sum;
61 add(sum, grad, -dsum*sum*sum, shape_x, grad);
62}
63
64void NURBS1DFiniteElement::CalcHessian (const IntegrationPoint &ip,
65 DenseMatrix &hessian) const
66{
67 Vector grad(dof);
68 Vector hess(hessian.Data(), dof);
69
70 kv[0]->CalcShape (shape_x, ijk[0], ip.x);
71 kv[0]->CalcDShape(grad, ijk[0], ip.x);
72 kv[0]->CalcD2Shape(hess, ijk[0], ip.x);
73
74 real_t sum = 0.0, dsum = 0.0, d2sum = 0.0;
75 for (int i = 0; i <= order; i++)
76 {
77 sum += (shape_x(i) *= weights(i));
78 dsum += ( grad(i) *= weights(i));
79 d2sum += ( hess(i) *= weights(i));
80 }
81
82 sum = 1.0/sum;
83 add(sum, hess, -2*dsum*sum*sum, grad, hess);
84 add(1.0, hess, (-d2sum + 2*dsum*dsum*sum)*sum*sum, shape_x, hess);
85}
86
87
88void NURBS2DFiniteElement::SetOrder() const
89{
90 orders[0] = kv[0]->GetOrder();
91 orders[1] = kv[1]->GetOrder();
92 shape_x.SetSize(orders[0]+1);
93 shape_y.SetSize(orders[1]+1);
94 dshape_x.SetSize(orders[0]+1);
95 dshape_y.SetSize(orders[1]+1);
96 d2shape_x.SetSize(orders[0]+1);
97 d2shape_y.SetSize(orders[1]+1);
98
99 order = max(orders[0], orders[1]);
100 dof = (orders[0] + 1)*(orders[1] + 1);
101 u.SetSize(dof);
102 du.SetSize(dof);
103 weights.SetSize(dof);
104}
105
106void NURBS2DFiniteElement::CalcShape(const IntegrationPoint &ip,
107 Vector &shape) const
108{
109 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
110 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
111
112 real_t sum = 0.0;
113 for (int o = 0, j = 0; j <= orders[1]; j++)
114 {
115 const real_t sy = shape_y(j);
116 for (int i = 0; i <= orders[0]; i++, o++)
117 {
118 sum += ( shape(o) = shape_x(i)*sy*weights(o) );
119 }
120 }
121
122 shape /= sum;
123}
124
125void NURBS2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
126 DenseMatrix &dshape) const
127{
128 real_t sum, dsum[2];
129
130 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
131 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
132
133 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
134 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
135
136 sum = dsum[0] = dsum[1] = 0.0;
137 for (int o = 0, j = 0; j <= orders[1]; j++)
138 {
139 const real_t sy = shape_y(j), dsy = dshape_y(j);
140 for (int i = 0; i <= orders[0]; i++, o++)
141 {
142 sum += ( u(o) = shape_x(i)*sy*weights(o) );
143
144 dsum[0] += ( dshape(o,0) = dshape_x(i)*sy *weights(o) );
145 dsum[1] += ( dshape(o,1) = shape_x(i)*dsy*weights(o) );
146 }
147 }
148
149 sum = 1.0/sum;
150 dsum[0] *= sum*sum;
151 dsum[1] *= sum*sum;
152
153 for (int o = 0; o < dof; o++)
154 {
155 dshape(o,0) = dshape(o,0)*sum - u(o)*dsum[0];
156 dshape(o,1) = dshape(o,1)*sum - u(o)*dsum[1];
157 }
158}
159
160void NURBS2DFiniteElement::CalcHessian (const IntegrationPoint &ip,
161 DenseMatrix &hessian) const
162{
163 real_t sum, dsum[2], d2sum[3];
164
165 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
166 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
167
168 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
169 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
170
171 kv[0]->CalcD2Shape(d2shape_x, ijk[0], ip.x);
172 kv[1]->CalcD2Shape(d2shape_y, ijk[1], ip.y);
173
174 sum = dsum[0] = dsum[1] = 0.0;
175 d2sum[0] = d2sum[1] = d2sum[2] = 0.0;
176 for (int o = 0, j = 0; j <= orders[1]; j++)
177 {
178 const real_t sy = shape_y(j), dsy = dshape_y(j), d2sy = d2shape_y(j);
179 for (int i = 0; i <= orders[0]; i++, o++)
180 {
181 const real_t sx = shape_x(i), dsx = dshape_x(i), d2sx = d2shape_x(i);
182 sum += ( u(o) = sx*sy*weights(o) );
183
184 dsum[0] += ( du(o,0) = dsx*sy*weights(o) );
185 dsum[1] += ( du(o,1) = sx*dsy*weights(o) );
186
187 d2sum[0] += ( hessian(o,0) = d2sx*sy*weights(o) );
188 d2sum[1] += ( hessian(o,1) = dsx*dsy*weights(o) );
189 d2sum[2] += ( hessian(o,2) = sx*d2sy*weights(o) );
190 }
191 }
192
193 sum = 1.0/sum;
194 dsum[0] *= sum;
195 dsum[1] *= sum;
196
197 d2sum[0] *= sum;
198 d2sum[1] *= sum;
199 d2sum[2] *= sum;
200
201 for (int o = 0; o < dof; o++)
202 {
203 hessian(o,0) = hessian(o,0)*sum
204 - 2*du(o,0)*sum*dsum[0]
205 + u[o]*sum*(2*dsum[0]*dsum[0] - d2sum[0]);
206
207 hessian(o,1) = hessian(o,1)*sum
208 - du(o,0)*sum*dsum[1]
209 - du(o,1)*sum*dsum[0]
210 + u[o]*sum*(2*dsum[0]*dsum[1] - d2sum[1]);
211
212 hessian(o,2) = hessian(o,2)*sum
213 - 2*du(o,1)*sum*dsum[1]
214 + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[2]);
215 }
216}
217
218
219void NURBS3DFiniteElement::SetOrder() const
220{
221 orders[0] = kv[0]->GetOrder();
222 orders[1] = kv[1]->GetOrder();
223 orders[2] = kv[2]->GetOrder();
224 shape_x.SetSize(orders[0]+1);
225 shape_y.SetSize(orders[1]+1);
226 shape_z.SetSize(orders[2]+1);
227
228 dshape_x.SetSize(orders[0]+1);
229 dshape_y.SetSize(orders[1]+1);
230 dshape_z.SetSize(orders[2]+1);
231
232 d2shape_x.SetSize(orders[0]+1);
233 d2shape_y.SetSize(orders[1]+1);
234 d2shape_z.SetSize(orders[2]+1);
235
236 order = max(max(orders[0], orders[1]), orders[2]);
237 dof = (orders[0] + 1)*(orders[1] + 1)*(orders[2] + 1);
238 u.SetSize(dof);
239 du.SetSize(dof);
240 weights.SetSize(dof);
241}
242
243void NURBS3DFiniteElement::CalcShape(const IntegrationPoint &ip,
244 Vector &shape) const
245{
246 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
247 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
248 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
249
250 real_t sum = 0.0;
251 for (int o = 0, k = 0; k <= orders[2]; k++)
252 {
253 const real_t sz = shape_z(k);
254 for (int j = 0; j <= orders[1]; j++)
255 {
256 const real_t sy_sz = shape_y(j)*sz;
257 for (int i = 0; i <= orders[0]; i++, o++)
258 {
259 sum += ( shape(o) = shape_x(i)*sy_sz*weights(o) );
260 }
261 }
262 }
263
264 shape /= sum;
265}
266
267void NURBS3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
268 DenseMatrix &dshape) const
269{
270 real_t sum, dsum[3];
271
272 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
273 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
274 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
275
276 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
277 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
278 kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
279
280 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
281 for (int o = 0, k = 0; k <= orders[2]; k++)
282 {
283 const real_t sz = shape_z(k), dsz = dshape_z(k);
284 for (int j = 0; j <= orders[1]; j++)
285 {
286 const real_t sy_sz = shape_y(j)* sz;
287 const real_t dsy_sz = dshape_y(j)* sz;
288 const real_t sy_dsz = shape_y(j)*dsz;
289 for (int i = 0; i <= orders[0]; i++, o++)
290 {
291 sum += ( u(o) = shape_x(i)*sy_sz*weights(o) );
292
293 dsum[0] += ( dshape(o,0) = dshape_x(i)* sy_sz *weights(o) );
294 dsum[1] += ( dshape(o,1) = shape_x(i)*dsy_sz *weights(o) );
295 dsum[2] += ( dshape(o,2) = shape_x(i)* sy_dsz*weights(o) );
296 }
297 }
298 }
299
300 sum = 1.0/sum;
301 dsum[0] *= sum*sum;
302 dsum[1] *= sum*sum;
303 dsum[2] *= sum*sum;
304
305 for (int o = 0; o < dof; o++)
306 {
307 dshape(o,0) = dshape(o,0)*sum - u(o)*dsum[0];
308 dshape(o,1) = dshape(o,1)*sum - u(o)*dsum[1];
309 dshape(o,2) = dshape(o,2)*sum - u(o)*dsum[2];
310 }
311}
312
313void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
314 DenseMatrix &hessian) const
315{
316 real_t sum, dsum[3], d2sum[6];
317
318 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
319 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
320 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
321
322 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
323 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
324 kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
325
326 kv[0]->CalcD2Shape(d2shape_x, ijk[0], ip.x);
327 kv[1]->CalcD2Shape(d2shape_y, ijk[1], ip.y);
328 kv[2]->CalcD2Shape(d2shape_z, ijk[2], ip.z);
329
330 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
331 d2sum[0] = d2sum[1] = d2sum[2] = d2sum[3] = d2sum[4] = d2sum[5] = 0.0;
332
333 for (int o = 0, k = 0; k <= orders[2]; k++)
334 {
335 const real_t sz = shape_z(k), dsz = dshape_z(k), d2sz = d2shape_z(k);
336 for (int j = 0; j <= orders[1]; j++)
337 {
338 const real_t sy = shape_y(j), dsy = dshape_y(j), d2sy = d2shape_y(j);
339 for (int i = 0; i <= orders[0]; i++, o++)
340 {
341 const real_t sx = shape_x(i), dsx = dshape_x(i), d2sx = d2shape_x(i);
342 sum += ( u(o) = sx*sy*sz*weights(o) );
343
344 dsum[0] += ( du(o,0) = dsx*sy*sz*weights(o) );
345 dsum[1] += ( du(o,1) = sx*dsy*sz*weights(o) );
346 dsum[2] += ( du(o,2) = sx*sy*dsz*weights(o) );
347
348 d2sum[0] += ( hessian(o,0) = d2sx*sy*sz*weights(o) );
349 d2sum[1] += ( hessian(o,1) = dsx*dsy*sz*weights(o) );
350 d2sum[2] += ( hessian(o,2) = dsx*sy*dsz*weights(o) );
351
352 d2sum[3] += ( hessian(o,3) = sx*dsy*dsz*weights(o) );
353
354 d2sum[4] += ( hessian(o,4) = sx*sy*d2sz*weights(o) );
355 d2sum[5] += ( hessian(o,5) = sx*d2sy*sz*weights(o) );
356 }
357 }
358 }
359
360 sum = 1.0/sum;
361 dsum[0] *= sum;
362 dsum[1] *= sum;
363 dsum[2] *= sum;
364
365 d2sum[0] *= sum;
366 d2sum[1] *= sum;
367 d2sum[2] *= sum;
368
369 d2sum[3] *= sum;
370 d2sum[4] *= sum;
371 d2sum[5] *= sum;
372
373 for (int o = 0; o < dof; o++)
374 {
375 hessian(o,0) = hessian(o,0)*sum
376 - 2*du(o,0)*sum*dsum[0]
377 + u[o]*sum*(2*dsum[0]*dsum[0] - d2sum[0]);
378
379 hessian(o,1) = hessian(o,1)*sum
380 - du(o,0)*sum*dsum[1]
381 - du(o,1)*sum*dsum[0]
382 + u[o]*sum*(2*dsum[0]*dsum[1] - d2sum[1]);
383
384 hessian(o,2) = hessian(o,2)*sum
385 - du(o,0)*sum*dsum[2]
386 - du(o,2)*sum*dsum[0]
387 + u[o]*sum*(2*dsum[0]*dsum[2] - d2sum[2]);
388
389 hessian(o,3) = hessian(o,3)*sum
390 - du(o,1)*sum*dsum[2]
391 - du(o,2)*sum*dsum[1]
392 + u[o]*sum*(2*dsum[1]*dsum[2] - d2sum[3]);
393
394 hessian(o,4) = hessian(o,4)*sum
395 - 2*du(o,2)*sum*dsum[2]
396 + u[o]*sum*(2*dsum[2]*dsum[2] - d2sum[4]);
397
398 hessian(o,5) = hessian(o,5)*sum
399 - 2*du(o,1)*sum*dsum[1]
400 + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[5]);
401
402 }
403}
404
405}
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseMatrix::Data
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:111
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
mfem::FiniteElement::dof
int dof
Number of degrees of freedom.
Definition fe_base.hpp:248
mfem::FiniteElement::orders
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition fe_base.hpp:250
mfem::FiniteElement::order
int order
Order/degree of the shape functions.
Definition fe_base.hpp:249
mfem::IntegrationPoint
Class for integration point with weight.
Definition intrules.hpp:35
mfem::NURBS1DFiniteElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_nurbs.cpp:31
mfem::NURBS1DFiniteElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_nurbs.cpp:45
mfem::NURBS1DFiniteElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_nurbs.cpp:64
mfem::NURBS1DFiniteElement::SetOrder
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition fe_nurbs.cpp:22
mfem::NURBS2DFiniteElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_nurbs.cpp:125
mfem::NURBS2DFiniteElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_nurbs.cpp:160
mfem::NURBS2DFiniteElement::SetOrder
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition fe_nurbs.cpp:88
mfem::NURBS2DFiniteElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_nurbs.cpp:106
mfem::NURBS3DFiniteElement::SetOrder
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition fe_nurbs.cpp:219
mfem::NURBS3DFiniteElement::CalcHessian
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_nurbs.cpp:313
mfem::NURBS3DFiniteElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_nurbs.cpp:243
mfem::NURBS3DFiniteElement::CalcDShape
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_nurbs.cpp:267
mfem::NURBSFiniteElement::kv
Array< const KnotVector * > kv
Definition fe_nurbs.hpp:26
mfem::Vector
Vector data type.
Definition vector.hpp:80
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
fe_nurbs.hpp
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
mfem::real_t
float real_t
Definition config.hpp:43
doxygen 1.11.0