MFEM  v4.6.0 Finite element discretization library
fe_nurbs.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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
17 namespace mfem
18 {
19
20 using namespace std;
21
23 {
24  order = kv[0]->GetOrder();
25  dof = order + 1;
26
27  weights.SetSize(dof);
28  shape_x.SetSize(dof);
29 }
30
32  Vector &shape) const
33 {
34  kv[0]->CalcShape(shape, ijk[0], ip.x);
35
36  double 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
46  DenseMatrix &dshape) const
47 {
49
50  kv[0]->CalcShape (shape_x, ijk[0], ip.x);
52
53  double 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;
62 }
63
65  DenseMatrix &hessian) const
66 {
68  Vector hess(hessian.Data(), dof);
69
70  kv[0]->CalcShape (shape_x, ijk[0], ip.x);
72  kv[0]->CalcD2Shape(hess, ijk[0], ip.x);
73
74  double 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;
84  add(1.0, hess, (-d2sum + 2*dsum*dsum*sum)*sum*sum, shape_x, hess);
85 }
86
87
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
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  double sum = 0.0;
113  for (int o = 0, j = 0; j <= orders[1]; j++)
114  {
115  const double 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
126  DenseMatrix &dshape) const
127 {
128  double 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 double 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
161  DenseMatrix &hessian) const
162 {
163  double 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 double 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 double 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
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
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  double sum = 0.0;
251  for (int o = 0, k = 0; k <= orders[2]; k++)
252  {
253  const double sz = shape_z(k);
254  for (int j = 0; j <= orders[1]; j++)
255  {
256  const double 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
268  DenseMatrix &dshape) const
269 {
270  double 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 double sz = shape_z(k), dsz = dshape_z(k);
284  for (int j = 0; j <= orders[1]; j++)
285  {
286  const double sy_sz = shape_y(j)* sz;
287  const double dsy_sz = dshape_y(j)* sz;
288  const double 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
314  DenseMatrix &hessian) const
315 {
316  double 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 double sz = shape_z(k), dsz = dshape_z(k), d2sz = d2shape_z(k);
336  for (int j = 0; j <= orders[1]; j++)
337  {
338  const double 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 double 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 }
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
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
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.cpp:219
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
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.cpp:22
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
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
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
STL namespace.
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
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
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
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.cpp:88
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
Class for integration point with weight.
Definition: intrules.hpp:31
Vector data type.
Definition: vector.hpp:58
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
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