MFEM  v4.5.2
Finite element discretization library
doftrans.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 "fem.hpp"
13 
14 namespace mfem
15 {
16 
18 {
20 }
21 
23 {
24  for (int c=0; c<V.Width(); c++)
25  {
27  }
28 }
29 
31 {
33 }
34 
36 {
39 }
40 
42 {
43  Vector row;
44  for (int r=0; r<V.Height(); r++)
45  {
46  V.GetRow(r, row);
47  TransformDual(row);
48  V.SetRow(r, row);
49  }
50 }
51 
53 {
54  for (int c=0; c<V.Width(); c++)
55  {
57  }
58 }
59 
61 {
63 }
64 
65 void TransformPrimal(const DofTransformation *ran_dof_trans,
66  const DofTransformation *dom_dof_trans,
67  DenseMatrix &elmat)
68 {
69  if (ran_dof_trans && dom_dof_trans)
70  {
71  ran_dof_trans->TransformPrimalCols(elmat);
72  dom_dof_trans->TransformDualRows(elmat);
73  }
74  else if (ran_dof_trans)
75  {
76  ran_dof_trans->TransformPrimalCols(elmat);
77  }
78  else if (dom_dof_trans)
79  {
80  dom_dof_trans->TransformDualRows(elmat);
81  }
82  else
83  {
84  // If both transformations are NULL this function should not be called
85  }
86 }
87 
89 {
91 }
92 
93 void TransformDual(const DofTransformation *ran_dof_trans,
94  const DofTransformation *dom_dof_trans,
95  DenseMatrix &elmat)
96 {
97  if (ran_dof_trans && dom_dof_trans)
98  {
99  ran_dof_trans->TransformDualCols(elmat);
100  dom_dof_trans->TransformDualRows(elmat);
101  }
102  else if (ran_dof_trans)
103  {
104  ran_dof_trans->TransformDualCols(elmat);
105  }
106  else if (dom_dof_trans)
107  {
108  dom_dof_trans->TransformDualRows(elmat);
109  }
110  else
111  {
112  // If both transformations are NULL this function should not be called
113  }
114 }
115 
117 {
118  int size = doftrans_->Size();
119 
120  if ((Ordering::Type)ordering_ == Ordering::byNODES || vdim_ == 1)
121  {
122  for (int i=0; i<vdim_; i++)
123  {
124  doftrans_->TransformPrimal(&v[i*size]);
125  }
126  }
127  else
128  {
129  Vector vec(size);
130  for (int i=0; i<vdim_; i++)
131  {
132  for (int j=0; j<size; j++)
133  {
134  vec(j) = v[j*vdim_+i];
135  }
136  doftrans_->TransformPrimal(vec);
137  for (int j=0; j<size; j++)
138  {
139  v[j*vdim_+i] = vec(j);
140  }
141  }
142  }
143 }
144 
146 {
147  int size = doftrans_->Height();
148 
149  if ((Ordering::Type)ordering_ == Ordering::byNODES)
150  {
151  for (int i=0; i<vdim_; i++)
152  {
153  doftrans_->InvTransformPrimal(&v[i*size]);
154  }
155  }
156  else
157  {
158  Vector vec(size);
159  for (int i=0; i<vdim_; i++)
160  {
161  for (int j=0; j<size; j++)
162  {
163  vec(j) = v[j*vdim_+i];
164  }
165  doftrans_->InvTransformPrimal(vec);
166  for (int j=0; j<size; j++)
167  {
168  v[j*vdim_+i] = vec(j);
169  }
170  }
171  }
172 }
173 
175 {
176  int size = doftrans_->Size();
177 
178  if ((Ordering::Type)ordering_ == Ordering::byNODES)
179  {
180  for (int i=0; i<vdim_; i++)
181  {
182  doftrans_->TransformDual(&v[i*size]);
183  }
184  }
185  else
186  {
187  Vector vec(size);
188  for (int i=0; i<vdim_; i++)
189  {
190  for (int j=0; j<size; j++)
191  {
192  vec(j) = v[j*vdim_+i];
193  }
194  doftrans_->TransformDual(vec);
195  for (int j=0; j<size; j++)
196  {
197  v[j*vdim_+i] = vec(j);
198  }
199  }
200  }
201 }
202 
204 {
205  int size = doftrans_->Size();
206 
207  if ((Ordering::Type)ordering_ == Ordering::byNODES)
208  {
209  for (int i=0; i<vdim_; i++)
210  {
211  doftrans_->InvTransformDual(&v[i*size]);
212  }
213  }
214  else
215  {
216  Vector vec(size);
217  for (int i=0; i<vdim_; i++)
218  {
219  for (int j=0; j<size; j++)
220  {
221  vec(j) = v[j*vdim_+i];
222  }
223  doftrans_->InvTransformDual(vec);
224  for (int j=0; j<size; j++)
225  {
226  v[j*vdim_+i] = vec(j);
227  }
228  }
229  }
230 }
231 
232 const double ND_DofTransformation::T_data[24] =
233 {
234  1.0, 0.0, 0.0, 1.0,
235  -1.0, -1.0, 0.0, 1.0,
236  0.0, 1.0, -1.0, -1.0,
237  1.0, 0.0, -1.0, -1.0,
238  -1.0, -1.0, 1.0, 0.0,
239  0.0, 1.0, 1.0, 0.0
240 };
241 
242 const DenseTensor ND_DofTransformation
243 ::T(const_cast<double*>(ND_DofTransformation::T_data), 2, 2, 6);
244 
245 const double ND_DofTransformation::TInv_data[24] =
246 {
247  1.0, 0.0, 0.0, 1.0,
248  -1.0, -1.0, 0.0, 1.0,
249  -1.0, -1.0, 1.0, 0.0,
250  1.0, 0.0, -1.0, -1.0,
251  0.0, 1.0, -1.0, -1.0,
252  0.0, 1.0, 1.0, 0.0
253 };
254 
255 const DenseTensor ND_DofTransformation
256 ::TInv(const_cast<double*>(TInv_data), 2, 2, 6);
257 
259  : DofTransformation(size)
260  , order(p)
261  , nedofs(p)
262  , nfdofs(p*(p-1))
263 {
264 }
265 
267  : ND_DofTransformation(p*(p + 2), p)
268 {
269 }
270 
272 {
273  // Return immediately when no face DoFs are present
274  if (nfdofs < 2) { return; }
275 
276  MFEM_VERIFY(Fo.Size() >= 1,
277  "Face orientations are unset in ND_TriDofTransformation");
278 
279  double data[2];
280  Vector v2(data, 2);
281 
282  // Transform face DoFs
283  for (int f=0; f<1; f++)
284  {
285  for (int i=0; i<nfdofs/2; i++)
286  {
287  v2 = &v[3*nedofs + f*nfdofs + 2*i];
288  T(Fo[f]).Mult(v2, &v[3*nedofs + f*nfdofs + 2*i]);
289  }
290  }
291 }
292 
293 void
295 {
296  // Return immediately when no face DoFs are present
297  if (nfdofs < 2) { return; }
298 
299  MFEM_VERIFY(Fo.Size() >= 1,
300  "Face orientations are unset in ND_TriDofTransformation");
301 
302  double data[2];
303  Vector v2(data, 2);
304 
305  // Transform face DoFs
306  for (int f=0; f<1; f++)
307  {
308  for (int i=0; i<nfdofs/2; i++)
309  {
310  v2 = &v[3*nedofs + f*nfdofs + 2*i];
311  TInv(Fo[f]).Mult(v2, &v[3*nedofs + f*nfdofs + 2*i]);
312  }
313  }
314 }
315 
316 void
318 {
319  // Return immediately when no face DoFs are present
320  if (nfdofs < 2) { return; }
321 
322  MFEM_VERIFY(Fo.Size() >= 1,
323  "Face orientations are unset in ND_TriDofTransformation");
324 
325  double data[2];
326  Vector v2(data, 2);
327 
328  // Transform face DoFs
329  for (int f=0; f<1; f++)
330  {
331  for (int i=0; i<nfdofs/2; i++)
332  {
333  v2 = &v[3*nedofs + f*nfdofs + 2*i];
334  TInv(Fo[f]).MultTranspose(v2, &v[3*nedofs + f*nfdofs + 2*i]);
335  }
336  }
337 }
338 
339 void
341 {
342  // Return immediately when no face DoFs are present
343  if (nfdofs < 2) { return; }
344 
345  MFEM_VERIFY(Fo.Size() >= 1,
346  "Face orientations are unset in ND_TriDofTransformation");
347 
348  double data[2];
349  Vector v2(data, 2);
350 
351  // Transform face DoFs
352  for (int f=0; f<1; f++)
353  {
354  for (int i=0; i<nfdofs/2; i++)
355  {
356  v2 = &v[3*nedofs + f*nfdofs + 2*i];
357  T(Fo[f]).MultTranspose(v2, &v[3*nedofs + f*nfdofs + 2*i]);
358  }
359  }
360 }
361 
363  : ND_DofTransformation(p*(p + 2)*(p + 3)/2, p)
364 {
365 }
366 
368 {
369  // Return immediately when no face DoFs are present
370  if (nfdofs < 2) { return; }
371 
372  MFEM_VERIFY(Fo.Size() >= 4,
373  "Face orientations are unset in ND_TetDofTransformation");
374 
375  double data[2];
376  Vector v2(data, 2);
377 
378  // Transform face DoFs
379  for (int f=0; f<4; f++)
380  {
381  for (int i=0; i<nfdofs/2; i++)
382  {
383  v2 = &v[6*nedofs + f*nfdofs + 2*i];
384  T(Fo[f]).Mult(v2, &v[6*nedofs + f*nfdofs + 2*i]);
385  }
386  }
387 }
388 
389 void
391 {
392  // Return immediately when no face DoFs are present
393  if (nfdofs < 2) { return; }
394 
395  MFEM_VERIFY(Fo.Size() >= 4,
396  "Face orientations are unset in ND_TetDofTransformation");
397 
398  double data[2];
399  Vector v2(data, 2);
400 
401  // Transform face DoFs
402  for (int f=0; f<4; f++)
403  {
404  for (int i=0; i<nfdofs/2; i++)
405  {
406  v2 = &v[6*nedofs + f*nfdofs + 2*i];
407  TInv(Fo[f]).Mult(v2, &v[6*nedofs + f*nfdofs + 2*i]);
408  }
409  }
410 }
411 
412 void
414 {
415  // Return immediately when no face DoFs are present
416  if (nfdofs < 2) { return; }
417 
418  MFEM_VERIFY(Fo.Size() >= 4,
419  "Face orientations are unset in ND_TetDofTransformation");
420 
421  double data[2];
422  Vector v2(data, 2);
423 
424  // Transform face DoFs
425  for (int f=0; f<4; f++)
426  {
427  for (int i=0; i<nfdofs/2; i++)
428  {
429  v2 = &v[6*nedofs + f*nfdofs + 2*i];
430  TInv(Fo[f]).MultTranspose(v2, &v[6*nedofs + f*nfdofs + 2*i]);
431  }
432  }
433 }
434 
435 void
437 {
438  // Return immediately when no face DoFs are present
439  if (nfdofs < 2) { return; }
440 
441  MFEM_VERIFY(Fo.Size() >= 4,
442  "Face orientations are unset in ND_TetDofTransformation");
443 
444  double data[2];
445  Vector v2(data, 2);
446 
447  // Transform face DoFs
448  for (int f=0; f<4; f++)
449  {
450  for (int i=0; i<nfdofs/2; i++)
451  {
452  v2 = &v[6*nedofs + f*nfdofs + 2*i];
453  T(Fo[f]).MultTranspose(v2, &v[6*nedofs + f*nfdofs + 2*i]);
454  }
455  }
456 }
457 
459  : ND_DofTransformation(3 * p * ((p + 1) * (p + 2))/2, p)
460 {
461 }
462 
464 {
465  // Return immediately when no face DoFs are present
466  if (nfdofs < 2) { return; }
467 
468  MFEM_VERIFY(Fo.Size() >= 2,
469  "Face orientations are unset in ND_WedgeDofTransformation");
470 
471  double data[2];
472  Vector v2(data, 2);
473 
474  // Transform triangular face DoFs
475  for (int f=0; f<2; f++)
476  {
477  for (int i=0; i<nfdofs/2; i++)
478  {
479  v2 = &v[9*nedofs + f*nfdofs + 2*i];
480  T(Fo[f]).Mult(v2, &v[9*nedofs + f*nfdofs + 2*i]);
481  }
482  }
483 }
484 
485 void
487 {
488  // Return immediately when no face DoFs are present
489  if (nfdofs < 2) { return; }
490 
491  MFEM_VERIFY(Fo.Size() >= 2,
492  "Face orientations are unset in ND_WedgeDofTransformation");
493 
494  double data[2];
495  Vector v2(data, 2);
496 
497  // Transform triangular face DoFs
498  for (int f=0; f<2; f++)
499  {
500  for (int i=0; i<nfdofs/2; i++)
501  {
502  v2 = &v[9*nedofs + f*nfdofs + 2*i];
503  TInv(Fo[f]).Mult(v2, &v[9*nedofs + f*nfdofs + 2*i]);
504  }
505  }
506 }
507 
508 void
510 {
511  // Return immediately when no face DoFs are present
512  if (nfdofs < 2) { return; }
513 
514  MFEM_VERIFY(Fo.Size() >= 2,
515  "Face orientations are unset in ND_WedgeDofTransformation");
516 
517  double data[2];
518  Vector v2(data, 2);
519 
520  // Transform triangular face DoFs
521  for (int f=0; f<2; f++)
522  {
523  for (int i=0; i<nfdofs/2; i++)
524  {
525  v2 = &v[9*nedofs + f*nfdofs + 2*i];
526  TInv(Fo[f]).MultTranspose(v2, &v[9*nedofs + f*nfdofs + 2*i]);
527  }
528  }
529 }
530 
531 void
533 {
534  // Return immediately when no face DoFs are present
535  if (nfdofs < 2) { return; }
536 
537  MFEM_VERIFY(Fo.Size() >= 2,
538  "Face orientations are unset in ND_WedgeDofTransformation");
539 
540  double data[2];
541  Vector v2(data, 2);
542 
543  // Transform triangular face DoFs
544  for (int f=0; f<2; f++)
545  {
546  for (int i=0; i<nfdofs/2; i++)
547  {
548  v2 = &v[9*nedofs + f*nfdofs + 2*i];
549  T(Fo[f]).MultTranspose(v2, &v[9*nedofs + f*nfdofs + 2*i]);
550  }
551  }
552 }
553 
554 } // namespace mfem
void SetRow(int r, const double *row)
Definition: densemat.cpp:2178
static const DenseTensor TInv
Definition: doftrans.hpp:215
void TransformPrimal(double *v) const
Definition: doftrans.cpp:271
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual void TransformPrimalCols(DenseMatrix &V) const
Transform groups of DoFs stored as dense matrices.
Definition: doftrans.cpp:22
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void InvTransformDual(double *v) const
Definition: doftrans.cpp:436
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:145
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:294
void InvTransformDual(double *v) const
Definition: doftrans.cpp:340
ND_DofTransformation(int size, int order)
Definition: doftrans.cpp:258
void InvTransformDual(double *v) const
Definition: doftrans.cpp:203
void TransformPrimal(double *v) const
Definition: doftrans.cpp:463
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void TransformPrimal(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition: doftrans.cpp:65
Type
Ordering methods:
Definition: fespace.hpp:33
static const double TInv_data[24]
Definition: doftrans.hpp:214
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1331
static const DenseTensor T
Definition: doftrans.hpp:215
void TransformDual(double *v) const
Definition: doftrans.cpp:174
void InvTransformDual(double *v) const
Definition: doftrans.cpp:532
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void TransformPrimal(double *v) const
Definition: doftrans.cpp:116
virtual void InvTransformPrimal(double *v) const =0
void TransformDual(double *v) const
Definition: doftrans.cpp:413
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:390
void TransformDual(double *v) const
Definition: doftrans.cpp:317
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1315
virtual void TransformPrimal(double *v) const =0
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:486
void TransformDual(double *v) const
Definition: doftrans.cpp:509
void TransformPrimal(double *v) const
Definition: doftrans.cpp:367
static const double T_data[24]
Definition: doftrans.hpp:213
virtual void TransformDualCols(DenseMatrix &V) const
Definition: doftrans.cpp:52
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
virtual void TransformDualRows(DenseMatrix &V) const
Transform groups of dual DoFs stored as dense matrices.
Definition: doftrans.cpp:41
Vector data type.
Definition: vector.hpp:60
void TransformDual(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition: doftrans.cpp:93
virtual void InvTransformDual(double *v) const =0
virtual void TransformDual(double *v) const =0