MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe_nurbs.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 // 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 
22 void 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 
31 void NURBS1DFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
45 void 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  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;
61  add(sum, grad, -dsum*sum*sum, shape_x, grad);
62 }
63 
64 void 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  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;
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 
88 void 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 
106 void 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  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 
125 void NURBS2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
160 void NURBS2DFiniteElement::CalcHessian (const IntegrationPoint &ip,
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 
219 void 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 
243 void 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  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 
267 void NURBS3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
313 void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
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 }
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::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::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::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
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::SetOrder
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.cpp:22
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.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
fe_nurbs.hpp
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::DenseMatrix::Data
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
mfem::NURBS3DFiniteElement::SetOrder
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.cpp:219
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
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::Vector
Vector data type.
Definition: vector.hpp:60
mfem::NURBS2DFiniteElement::SetOrder
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.cpp:88
mfem::u
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
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::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
Generated on Wed Dec 7 2022 10:41:45 for MFEM by   doxygen 1.8.5