MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
doftrans.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 #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  int nedofs = order; // number of DoFs per edge
343  int nfdofs = order*(order-1); // number of DoFs per face
344 
345  double data[2];
346  Vector v2(data, 2);
347 
348  // Transform face DoFs
349  for (int f=0; f<1; f++)
350  {
351  for (int i=0; i<nfdofs/2; i++)
352  {
353  v2 = &v[3*nedofs + f*nfdofs + 2*i];
354  T(Fo[f]).MultTranspose(v2, &v[3*nedofs + f*nfdofs + 2*i]);
355  }
356  }
357 }
358 
360  : ND_DofTransformation(p*(p + 2)*(p + 3)/2, p)
361 {
362 }
363 
365 {
366  // Return immediately when no face DoFs are present
367  if (nfdofs < 2) { return; }
368 
369  MFEM_VERIFY(Fo.Size() >= 4,
370  "Face orientations are unset in ND_TetDofTransformation");
371 
372  double data[2];
373  Vector v2(data, 2);
374 
375  // Transform face DoFs
376  for (int f=0; f<4; f++)
377  {
378  for (int i=0; i<nfdofs/2; i++)
379  {
380  v2 = &v[6*nedofs + f*nfdofs + 2*i];
381  T(Fo[f]).Mult(v2, &v[6*nedofs + f*nfdofs + 2*i]);
382  }
383  }
384 }
385 
386 void
388 {
389  // Return immediately when no face DoFs are present
390  if (nfdofs < 2) { return; }
391 
392  MFEM_VERIFY(Fo.Size() >= 4,
393  "Face orientations are unset in ND_TetDofTransformation");
394 
395  double data[2];
396  Vector v2(data, 2);
397 
398  // Transform face DoFs
399  for (int f=0; f<4; f++)
400  {
401  for (int i=0; i<nfdofs/2; i++)
402  {
403  v2 = &v[6*nedofs + f*nfdofs + 2*i];
404  TInv(Fo[f]).Mult(v2, &v[6*nedofs + f*nfdofs + 2*i]);
405  }
406  }
407 }
408 
409 void
411 {
412  // Return immediately when no face DoFs are present
413  if (nfdofs < 2) { return; }
414 
415  MFEM_VERIFY(Fo.Size() >= 4,
416  "Face orientations are unset in ND_TetDofTransformation");
417 
418  double data[2];
419  Vector v2(data, 2);
420 
421  // Transform face DoFs
422  for (int f=0; f<4; f++)
423  {
424  for (int i=0; i<nfdofs/2; i++)
425  {
426  v2 = &v[6*nedofs + f*nfdofs + 2*i];
427  TInv(Fo[f]).MultTranspose(v2, &v[6*nedofs + f*nfdofs + 2*i]);
428  }
429  }
430 }
431 
432 void
434 {
435  int nedofs = order; // number of DoFs per edge
436  int nfdofs = order*(order-1); // number of DoFs per face
437 
438  double data[2];
439  Vector v2(data, 2);
440 
441  // Transform face DoFs
442  for (int f=0; f<4; f++)
443  {
444  for (int i=0; i<nfdofs/2; i++)
445  {
446  v2 = &v[6*nedofs + f*nfdofs + 2*i];
447  T(Fo[f]).MultTranspose(v2, &v[6*nedofs + f*nfdofs + 2*i]);
448  }
449  }
450 }
451 
453  : ND_DofTransformation(3 * p * ((p + 1) * (p + 2))/2, p)
454 {
455 }
456 
458 {
459  // Return immediately when no face DoFs are present
460  if (nfdofs < 2) { return; }
461 
462  MFEM_VERIFY(Fo.Size() >= 2,
463  "Face orientations are unset in ND_WedgeDofTransformation");
464 
465  double data[2];
466  Vector v2(data, 2);
467 
468  // Transform triangular face DoFs
469  for (int f=0; f<2; f++)
470  {
471  for (int i=0; i<nfdofs/2; i++)
472  {
473  v2 = &v[9*nedofs + f*nfdofs + 2*i];
474  T(Fo[f]).Mult(v2, &v[9*nedofs + f*nfdofs + 2*i]);
475  }
476  }
477 }
478 
479 void
481 {
482  // Return immediately when no face DoFs are present
483  if (nfdofs < 2) { return; }
484 
485  MFEM_VERIFY(Fo.Size() >= 2,
486  "Face orientations are unset in ND_WedgeDofTransformation");
487 
488  double data[2];
489  Vector v2(data, 2);
490 
491  // Transform triangular face DoFs
492  for (int f=0; f<2; f++)
493  {
494  for (int i=0; i<nfdofs/2; i++)
495  {
496  v2 = &v[9*nedofs + f*nfdofs + 2*i];
497  TInv(Fo[f]).Mult(v2, &v[9*nedofs + f*nfdofs + 2*i]);
498  }
499  }
500 }
501 
502 void
504 {
505  // Return immediately when no face DoFs are present
506  if (nfdofs < 2) { return; }
507 
508  MFEM_VERIFY(Fo.Size() >= 2,
509  "Face orientations are unset in ND_WedgeDofTransformation");
510 
511  double data[2];
512  Vector v2(data, 2);
513 
514  // Transform triangular face DoFs
515  for (int f=0; f<2; f++)
516  {
517  for (int i=0; i<nfdofs/2; i++)
518  {
519  v2 = &v[9*nedofs + f*nfdofs + 2*i];
520  TInv(Fo[f]).MultTranspose(v2, &v[9*nedofs + f*nfdofs + 2*i]);
521  }
522  }
523 }
524 
525 void
527 {
528  // Return immediately when no face DoFs are present
529  if (nfdofs < 2) { return; }
530 
531  MFEM_VERIFY(Fo.Size() >= 2,
532  "Face orientations are unset in ND_WedgeDofTransformation");
533 
534  double data[2];
535  Vector v2(data, 2);
536 
537  // Transform triangular face DoFs
538  for (int f=0; f<2; f++)
539  {
540  for (int i=0; i<nfdofs/2; i++)
541  {
542  v2 = &v[9*nedofs + f*nfdofs + 2*i];
543  T(Fo[f]).MultTranspose(v2, &v[9*nedofs + f*nfdofs + 2*i]);
544  }
545  }
546 }
547 
548 } // namespace mfem
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:294
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void TransformDual(double *v) const
Definition: doftrans.cpp:410
void SetRow(int r, const double *row)
Definition: densemat.cpp:1773
static const DenseTensor TInv
Definition: doftrans.hpp:215
void TransformPrimal(double *v) const
Definition: doftrans.cpp:457
void InvTransformDual(double *v) const
Definition: doftrans.cpp:526
int Height() const
Definition: doftrans.hpp:69
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void TransformDual(double *v) const
Definition: doftrans.cpp:174
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
ND_DofTransformation(int size, int order)
Definition: doftrans.cpp:258
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void InvTransformDual(double *v) const
Definition: doftrans.cpp:433
void InvTransformDual(double *v) const
Definition: doftrans.cpp:203
virtual void TransformPrimalCols(DenseMatrix &V) const
Transform groups of DoFs stored as dense matrices.
Definition: doftrans.cpp:22
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void TransformDual(double *v) const
Definition: doftrans.cpp:503
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void TransformPrimal(const DofTransformation *ran_dof_trans, const DofTransformation *dom_dof_trans, DenseMatrix &elmat)
Definition: doftrans.cpp:65
void TransformDual(double *v) const
Definition: doftrans.cpp:317
Type
Ordering methods:
Definition: fespace.hpp:33
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1254
virtual void TransformDualCols(DenseMatrix &V) const
Definition: doftrans.cpp:52
static const double TInv_data[24]
Definition: doftrans.hpp:214
void TransformPrimal(double *v) const
Definition: doftrans.cpp:364
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1270
static const DenseTensor T
Definition: doftrans.hpp:215
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:145
virtual void InvTransformPrimal(double *v) const =0
virtual void TransformPrimal(double *v) const =0
void TransformPrimal(double *v) const
Definition: doftrans.cpp:271
virtual void TransformDualRows(DenseMatrix &V) const
Transform groups of dual DoFs stored as dense matrices.
Definition: doftrans.cpp:41
static const double T_data[24]
Definition: doftrans.hpp:213
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
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:480
void InvTransformPrimal(double *v) const
Definition: doftrans.cpp:387
void InvTransformDual(double *v) const
Definition: doftrans.cpp:340
virtual void InvTransformDual(double *v) const =0
virtual void TransformDual(double *v) const =0
void TransformPrimal(double *v) const
Definition: doftrans.cpp:116