Compare commits
7
Commits
MA
...
blf-full-op-dev
| Author | SHA1 | Date | |
|---|---|---|---|
 Stowell, Mark L
|
de20f03613 | ||
|
Stowell, Mark L
|
d58d632796 | ||
|
Stowell, Mark L
|
a9f8e29727 | ||
|
Stowell, Mark L
|
9a78d866cf | ||
|
Stowell, Mark L
|
51b302336c | ||
|
Stowell, Mark L
|
faf53d37f9 | ||
|
Stowell, Mark L
|
fc561569d1 |
@@ -1766,8 +1766,13 @@ void MixedBilinearForm::FormRectangularSystemMatrix(
|
||||
mat_e = new SparseMatrix(mat->Height(), mat->Width());
|
||||
mat->EliminateCols(ess_trial_tdof_marker, *mat_e);
|
||||
|
||||
Array<int> cols;
|
||||
Vector srow;
|
||||
|
||||
for (int i=0; i<test_tdof_list.Size(); ++i)
|
||||
{
|
||||
mat->GetRow(test_tdof_list[i], cols, srow);
|
||||
mat_e->AddRow(test_tdof_list[i], cols, srow);
|
||||
mat->EliminateRow(test_tdof_list[i]);
|
||||
}
|
||||
mat_e->Finalize();
|
||||
|
||||
+68
-21
@@ -276,9 +276,19 @@ public:
|
||||
|
||||
/** @brief Add the original uneliminated matrix vector multiple to a vector.
|
||||
The original matrix is \f$ M + Me \f$ so we have:
|
||||
\f$ y += M x + M_e x \f$ */
|
||||
void FullAddMult(const Vector &x, Vector &y) const
|
||||
{ mat->AddMult(x, y); mat_e->AddMult(x, y); }
|
||||
\f$ y += a * M x + a * M_e x \f$ */
|
||||
void FullAddMult(const Vector &x, Vector &y, const double a = 1.0) const
|
||||
{ mat->AddMult(x, y, a); mat_e->AddMult(x, y, a); }
|
||||
|
||||
/// Matrix transpose vector multiplication: \f$ y = M^T x \f$
|
||||
virtual void MultTranspose(const Vector & x, Vector & y) const
|
||||
{ y = 0.0; AddMultTranspose (x, y); }
|
||||
|
||||
/** @brief Matrix transpose vector multiplication with the original
|
||||
uneliminated matrix. The original matrix is \f$ M + M_e \f$ so we have:
|
||||
\f$ y = M^T x + {M_e}^T x \f$ */
|
||||
void FullMultTranspose(const Vector &x, Vector &y) const
|
||||
{ mat->MultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
|
||||
|
||||
/// Add the matrix transpose vector multiplication: \f$ y += a M^T x \f$
|
||||
virtual void AddMultTranspose(const Vector & x, Vector & y,
|
||||
@@ -287,18 +297,19 @@ public:
|
||||
|
||||
/** @brief Add the original uneliminated matrix transpose vector
|
||||
multiple to a vector. The original matrix is \f$ M + M_e \f$
|
||||
so we have: \f$ y += M^T x + {M_e}^T x \f$ */
|
||||
void FullAddMultTranspose(const Vector & x, Vector & y) const
|
||||
{ mat->AddMultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
|
||||
|
||||
/// Matrix transpose vector multiplication: \f$ y = M^T x \f$
|
||||
virtual void MultTranspose(const Vector & x, Vector & y) const
|
||||
{ y = 0.0; AddMultTranspose (x, y); }
|
||||
so we have: \f$ y += a * M^T x + a * {M_e}^T x \f$ */
|
||||
void FullAddMultTranspose(const Vector & x, Vector & y,
|
||||
const double a = 1.0) const
|
||||
{ mat->AddMultTranspose(x, y, a); mat_e->AddMultTranspose(x, y, a); }
|
||||
|
||||
/// Compute \f$ y^T M x \f$
|
||||
double InnerProduct(const Vector &x, const Vector &y) const
|
||||
{ return mat->InnerProduct (x, y); }
|
||||
|
||||
/// Compute inner product for full uneliminated matrix \f$ y^T M x + y^T M_e x \f$
|
||||
double FullInnerProduct(const Vector &x, const Vector &y) const
|
||||
{ return mat->InnerProduct(x, y) + mat_e->InnerProduct(x, y); }
|
||||
|
||||
/// Returns a pointer to (approximation) of the matrix inverse: \f$ M^{-1} \f$
|
||||
virtual MatrixInverse *Inverse() const;
|
||||
|
||||
@@ -434,8 +445,14 @@ public:
|
||||
recovered by calling RecoverFEMSolution() (with the same vectors @a X,
|
||||
@a b, and @a x).
|
||||
|
||||
NOTE: If there are no transformations, @a X simply reuses the data of
|
||||
@a x. */
|
||||
@note If there are no transformations, @a X simply reuses the data of
|
||||
@a x.
|
||||
|
||||
@note This method does modify the bilinear form operator. For example,
|
||||
calls to Mult() will produce different results before and after
|
||||
use of this method. Use FullMult() to obtain the original behavior.
|
||||
Similar methods exist for AddMult(), MultTranspose(), etc..
|
||||
*/
|
||||
virtual void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
|
||||
Vector &b, OperatorHandle &A, Vector &X,
|
||||
Vector &B, int copy_interior = 0);
|
||||
@@ -590,10 +607,6 @@ public:
|
||||
void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x,
|
||||
Vector &b);
|
||||
|
||||
/// Compute inner product for full uneliminated matrix \f$ y^T M x + y^T M_e x \f$
|
||||
double FullInnerProduct(const Vector &x, const Vector &y) const
|
||||
{ return mat->InnerProduct(x, y) + mat_e->InnerProduct(x, y); }
|
||||
|
||||
/// Update the @a FiniteElementSpace and delete all data associated with the old one.
|
||||
virtual void Update(FiniteElementSpace *nfes = NULL);
|
||||
|
||||
@@ -714,13 +727,42 @@ public:
|
||||
/// Matrix multiplication: \f$ y = M x \f$
|
||||
virtual void Mult(const Vector & x, Vector & y) const;
|
||||
|
||||
/** @brief Matrix vector multiplication with the original uneliminated
|
||||
matrix. The original matrix is \f$ M + M_e \f$ so we have:
|
||||
\f$ y = M x + M_e x \f$ */
|
||||
void FullMult(const Vector &x, Vector &y) const
|
||||
{ mat->Mult(x, y); mat_e->AddMult(x, y); }
|
||||
|
||||
/// Add the matrix vector multiple to a vector: \f$ y += a M x \f$
|
||||
virtual void AddMult(const Vector & x, Vector & y,
|
||||
const double a = 1.0) const;
|
||||
|
||||
/** @brief Add the original uneliminated matrix vector multiple to a vector.
|
||||
The original matrix is \f$ M + Me \f$ so we have:
|
||||
\f$ y += a * M x + a * M_e x \f$ */
|
||||
void FullAddMult(const Vector &x, Vector &y, const double a = 1.0) const
|
||||
{ mat->AddMult(x, y, a); mat_e->AddMult(x, y, a); }
|
||||
|
||||
/// Matrix transpose vector multiplication: \f$ y = M^T x \f$
|
||||
virtual void MultTranspose(const Vector & x, Vector & y) const;
|
||||
|
||||
/** @brief Matrix transpose vector multiplication with the original
|
||||
uneliminated matrix. The original matrix is \f$ M + M_e \f$ so we have:
|
||||
\f$ y = M^T x + {M_e}^T x \f$ */
|
||||
void FullMultTranspose(const Vector &x, Vector &y) const
|
||||
{ mat->MultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
|
||||
|
||||
/// Add the matrix transpose vector multiplication: \f$ y += a M^T x \f$
|
||||
virtual void AddMultTranspose(const Vector & x, Vector & y,
|
||||
const double a = 1.0) const;
|
||||
|
||||
/** @brief Add the original uneliminated matrix transpose vector
|
||||
multiple to a vector. The original matrix is \f$ M + M_e \f$
|
||||
so we have: \f$ y += a * M^T x + a * {M_e}^T x \f$ */
|
||||
void FullAddMultTranspose(const Vector & x, Vector & y,
|
||||
const double a = 1.0) const
|
||||
{ mat->AddMultTranspose(x, y, a); mat_e->AddMultTranspose(x, y, a); }
|
||||
|
||||
virtual MatrixInverse *Inverse() const;
|
||||
|
||||
/// Finalizes the matrix initialization.
|
||||
@@ -904,12 +946,17 @@ public:
|
||||
A.MakeRef(*A_ptr);
|
||||
}
|
||||
|
||||
/** @brief Form the linear system A X = B, corresponding to this mixed bilinear
|
||||
form and the linear form @a b(.).
|
||||
/** @brief Form the linear system A X = B, corresponding to this mixed
|
||||
bilinear form and the linear form @a b(.). */
|
||||
/** Return in @a A a *reference* to the system matrix that is
|
||||
column-constrained. The reference will be invalidated when
|
||||
SetOperatorType(), Update(), or the destructor is called.
|
||||
|
||||
Return in @a A a *reference* to the system matrix that is column-constrained.
|
||||
The reference will be invalidated when SetOperatorType(), Update(), or the
|
||||
destructor is called. */
|
||||
@note This method does modify the bilinear form operator. For example,
|
||||
calls to Mult() will produce different results before and after
|
||||
use of this method. Use FullMult() to obtain the original behavior.
|
||||
Similar methods exist for AddMult(), MultTranspose(), etc..
|
||||
*/
|
||||
virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
|
||||
const Array<int> &test_tdof_list,
|
||||
Vector &x, Vector &b,
|
||||
|
||||
@@ -47,6 +47,7 @@ set(UNIT_TESTS_SRCS
|
||||
mesh/test_pmesh.cpp
|
||||
mesh/test_periodic_mesh.cpp
|
||||
mesh/test_vtu.cpp
|
||||
fem/common_get_mesh.cpp
|
||||
fem/test_1d_bilininteg.cpp
|
||||
fem/test_2d_bilininteg.cpp
|
||||
fem/test_3d_bilininteg.cpp
|
||||
|
||||
@@ -0,0 +1,416 @@
|
||||
// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
|
||||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||||
//
|
||||
// This file is part of the MFEM library. For more information and source code
|
||||
// availability visit https://mfem.org.
|
||||
//
|
||||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||||
// CONTRIBUTING.md for details.
|
||||
|
||||
#include "common_get_mesh.hpp"
|
||||
|
||||
using namespace mfem;
|
||||
|
||||
namespace mfem_test_fem
|
||||
{
|
||||
|
||||
Mesh * GetMesh(MeshType type, double lx, double ly, double lz)
|
||||
{
|
||||
Mesh * mesh = NULL;
|
||||
double c[3];
|
||||
int v[8];
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case SEGMENT:
|
||||
mesh = new Mesh(1, 2, 1);
|
||||
c[0] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx;
|
||||
mesh->AddVertex(c);
|
||||
v[0] = 0; v[1] = 1;
|
||||
mesh->AddSegment(v);
|
||||
{
|
||||
Element * el = mesh->NewElement(Geometry::POINT);
|
||||
el->SetAttribute(1);
|
||||
el->SetVertices(&v[0]);
|
||||
mesh->AddBdrElement(el);
|
||||
}
|
||||
{
|
||||
Element * el = mesh->NewElement(Geometry::POINT);
|
||||
el->SetAttribute(2);
|
||||
el->SetVertices(&v[1]);
|
||||
mesh->AddBdrElement(el);
|
||||
}
|
||||
break;
|
||||
case QUADRILATERAL:
|
||||
mesh = new Mesh(2, 4, 1);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
|
||||
mesh->AddQuad(v);
|
||||
break;
|
||||
case TRIANGLE2A:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 0;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case TRIANGLE2B:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 1; v[1] = 2; v[2] = 0;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 2;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case TRIANGLE2C:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 2; v[1] = 0; v[2] = 1;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 0; v[1] = 2; v[2] = 3;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case TRIANGLE4:
|
||||
mesh = new Mesh(2, 5, 4);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = 0.5 * ly;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 1; v[1] = 2; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case MIXED2D:
|
||||
mesh = new Mesh(2, 6, 4);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * ly; c[1] = 0.5 * ly;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx - 0.5 * ly; c[1] = 0.5 * ly;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 5; v[3] = 4;
|
||||
mesh->AddQuad(v);
|
||||
v[0] = 1; v[1] = 2; v[2] = 5;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 5;
|
||||
mesh->AddQuad(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case HEXAHEDRON:
|
||||
mesh = new Mesh(3, 8, 1);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
|
||||
v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2A:
|
||||
case HEXAHEDRON2B:
|
||||
case HEXAHEDRON2C:
|
||||
case HEXAHEDRON2D:
|
||||
mesh = new Mesh(3, 12, 2);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 5; v[2] = 11; v[3] = 6;
|
||||
v[4] = 1; v[5] = 4; v[6] = 10; v[7] = 7;
|
||||
mesh->AddHex(v);
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case HEXAHEDRON2A: // Face Orientation 1
|
||||
v[0] = 4; v[1] = 10; v[2] = 7; v[3] = 1;
|
||||
v[4] = 3; v[5] = 9; v[6] = 8; v[7] = 2;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2B: // Face Orientation 3
|
||||
v[0] = 10; v[1] = 7; v[2] = 1; v[3] = 4;
|
||||
v[4] = 9; v[5] = 8; v[6] = 2; v[7] = 3;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2C: // Face Orientation 5
|
||||
v[0] = 7; v[1] = 1; v[2] = 4; v[3] = 10;
|
||||
v[4] = 8; v[5] = 2; v[6] = 3; v[7] = 9;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2D: // Face Orientation 7
|
||||
v[0] = 1; v[1] = 4; v[2] = 10; v[3] = 7;
|
||||
v[4] = 2; v[5] = 3; v[6] = 9; v[7] = 8;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
default:
|
||||
// Cannot happen
|
||||
break;
|
||||
}
|
||||
break;
|
||||
case WEDGE2:
|
||||
mesh = new Mesh(3, 8, 2);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 4; v[4] = 5; v[5] = 6;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 0; v[1] = 2; v[2] = 3; v[3] = 4; v[4] = 6; v[5] = 7;
|
||||
mesh->AddWedge(v);
|
||||
break;
|
||||
case TETRAHEDRA:
|
||||
mesh = new Mesh(3, 8, 5);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 2; v[2] = 7; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 6; v[1] = 7; v[2] = 2; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 4; v[1] = 7; v[2] = 5; v[3] = 0;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 1; v[1] = 0; v[2] = 5; v[3] = 2;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 3; v[1] = 7; v[2] = 0; v[3] = 2;
|
||||
mesh->AddTet(v);
|
||||
break;
|
||||
case WEDGE4:
|
||||
mesh = new Mesh(3, 10, 4);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = 0.5 * ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lx; c[1] = 0.5 * ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 4; v[3] = 5; v[4] = 6; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 1; v[1] = 2; v[2] = 4; v[3] = 6; v[4] = 7; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 7; v[4] = 8; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 4; v[3] = 8; v[4] = 5; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
break;
|
||||
case MIXED3D6:
|
||||
mesh = new Mesh(3, 12, 6);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lz; c[1] = 0.5 * lz; c[2] = 0.5 * lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx - 0.5 * lz; c[1] = 0.5 * lz; c[2] = 0.5 * lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx - 0.5 * lz; c[1] = ly - 0.5 * lz; c[2] = 0.5 * lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * lz; c[1] = ly - 0.5 * lz; c[2] = 0.5 * lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
|
||||
v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
|
||||
mesh->AddHex(v);
|
||||
v[0] = 0; v[1] = 4; v[2] = 8; v[3] = 1; v[4] = 5; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 1; v[1] = 5; v[2] = 9; v[3] = 2; v[4] = 6; v[5] = 10;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 2; v[1] = 6; v[2] = 10; v[3] = 3; v[4] = 7; v[5] = 11;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 3; v[1] = 7; v[2] = 11; v[3] = 0; v[4] = 4; v[5] = 8;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 4; v[1] = 5; v[2] = 6; v[3] = 7;
|
||||
v[4] = 8; v[5] = 9; v[6] = 10; v[7] = 11;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case MIXED3D8:
|
||||
mesh = new Mesh(3, 10, 8);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
c[0] = 0.25 * lx; c[1] = 0.5 * ly; c[2] = 0.5 * lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.75 * lx; c[1] = 0.5 * ly; c[2] = 0.5 * lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = 0.0; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = lx; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = ly; c[2] = lz;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 3; v[2] = 4; v[3] = 1; v[4] = 2; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 3; v[1] = 9; v[2] = 4; v[3] = 2; v[4] = 8; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 9; v[1] = 6; v[2] = 4; v[3] = 8; v[4] = 7; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 6; v[1] = 0; v[2] = 4; v[3] = 7; v[4] = 1; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 0; v[1] = 3; v[2] = 9; v[3] = 4;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 0; v[1] = 9; v[2] = 6; v[3] = 4;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 1; v[1] = 7; v[2] = 2; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 8; v[1] = 2; v[2] = 7; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
break;
|
||||
}
|
||||
mesh->FinalizeTopology();
|
||||
|
||||
return mesh;
|
||||
}
|
||||
|
||||
} // namespace mfem_test_fem
|
||||
@@ -0,0 +1,41 @@
|
||||
// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
|
||||
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||||
// LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||||
//
|
||||
// This file is part of the MFEM library. For more information and source code
|
||||
// availability visit https://mfem.org.
|
||||
//
|
||||
// MFEM is free software; you can redistribute it and/or modify it under the
|
||||
// terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||||
// CONTRIBUTING.md for details.
|
||||
|
||||
#include "mfem.hpp"
|
||||
|
||||
namespace mfem_test_fem
|
||||
{
|
||||
|
||||
enum MeshType
|
||||
{
|
||||
SEGMENT = 0,
|
||||
QUADRILATERAL = 1,
|
||||
TRIANGLE2A = 2,
|
||||
TRIANGLE2B = 3,
|
||||
TRIANGLE2C = 4,
|
||||
TRIANGLE4 = 5,
|
||||
MIXED2D = 6,
|
||||
HEXAHEDRON = 7,
|
||||
HEXAHEDRON2A = 8,
|
||||
HEXAHEDRON2B = 9,
|
||||
HEXAHEDRON2C = 10,
|
||||
HEXAHEDRON2D = 11,
|
||||
WEDGE2 = 12,
|
||||
TETRAHEDRA = 13,
|
||||
WEDGE4 = 14,
|
||||
MIXED3D6 = 15,
|
||||
MIXED3D8 = 16
|
||||
};
|
||||
|
||||
mfem::Mesh * GetMesh(MeshType type,
|
||||
double lx = 1.0, double ly = 1.0, double lz = 1.0);
|
||||
|
||||
}
|
||||
@@ -11,10 +11,19 @@
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "unit_tests.hpp"
|
||||
#include "common_get_mesh.hpp"
|
||||
|
||||
#include <iostream>
|
||||
|
||||
using namespace mfem;
|
||||
using namespace mfem_test_fem;
|
||||
|
||||
namespace bilinearform
|
||||
{
|
||||
|
||||
static double a_ = 5.0;
|
||||
static double b_ = 3.0;
|
||||
static double c_ = 2.0;
|
||||
|
||||
TEST_CASE("Test order of boundary integrators",
|
||||
"[BilinearForm]")
|
||||
@@ -142,3 +151,293 @@ TEST_CASE("FormLinearSystem/SolutionScope",
|
||||
REQUIRE(AsConst(sol)(bdr_dof) == 0.0);
|
||||
}
|
||||
}
|
||||
|
||||
enum FEType
|
||||
{
|
||||
H1_FEC = 0,
|
||||
ND_FEC,
|
||||
RT_FEC,
|
||||
L2V_FEC,
|
||||
L2I_FEC,
|
||||
};
|
||||
|
||||
TEST_CASE("BilinearForm Full Ops",
|
||||
"[BilinearForm]")
|
||||
{
|
||||
int order = 2;
|
||||
double alpha = M_E;
|
||||
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
mesh->UniformRefinement();
|
||||
|
||||
Vector oneVec(dim); oneVec = 1.0;
|
||||
|
||||
ConstantCoefficient oneCoef(1.0);
|
||||
VectorConstantCoefficient oneVecCoef(oneVec);
|
||||
|
||||
for (int ft = (int)FEType::H1_FEC; ft <= (int)FEType::RT_FEC; ft++)
|
||||
{
|
||||
// if (ft == (int)FEType::ND_FEC || ft == (int)FEType::RT_FEC)
|
||||
// { continue; }
|
||||
bool vec = (ft == (int)FEType::ND_FEC || ft == (int)FEType::RT_FEC);
|
||||
|
||||
if (dim == 1 && vec) { continue; }
|
||||
if (vec && (mt == (int)MeshType::WEDGE2 ||
|
||||
mt == (int)MeshType::WEDGE4 ||
|
||||
mt == (int)MeshType::MIXED3D6 ||
|
||||
mt == (int)MeshType::MIXED3D8))
|
||||
{ continue; }
|
||||
|
||||
SECTION("Integral of field " + std::to_string(ft) +
|
||||
" on mesh type " + std::to_string(mt) )
|
||||
{
|
||||
|
||||
FiniteElementCollection *fec = NULL;
|
||||
switch ((FEType)ft)
|
||||
{
|
||||
case FEType::H1_FEC:
|
||||
fec = new H1_FECollection(order, dim);
|
||||
break;
|
||||
case FEType::ND_FEC:
|
||||
fec = new ND_FECollection(order, dim);
|
||||
break;
|
||||
case FEType::RT_FEC:
|
||||
fec = new RT_FECollection(order-1, dim);
|
||||
break;
|
||||
case FEType::L2V_FEC:
|
||||
fec = new L2_FECollection(order-1, dim);
|
||||
break;
|
||||
case FEType::L2I_FEC:
|
||||
fec = new L2_FECollection(order, dim,
|
||||
BasisType::GaussLegendre,
|
||||
FiniteElement::INTEGRAL);
|
||||
break;
|
||||
default:
|
||||
MFEM_ABORT("Invalid vector FE type");
|
||||
}
|
||||
FiniteElementSpace fespace(mesh, fec);
|
||||
|
||||
Array<int> ess_tdof_list;
|
||||
if (mesh->bdr_attributes.Size())
|
||||
{
|
||||
Array<int> ess_bdr(mesh->bdr_attributes.Max());
|
||||
ess_bdr = 1;
|
||||
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
|
||||
}
|
||||
|
||||
GridFunction u(&fespace);
|
||||
if (!vec)
|
||||
{
|
||||
u.ProjectCoefficient(oneCoef);
|
||||
}
|
||||
else
|
||||
{
|
||||
u.ProjectCoefficient(oneVecCoef);
|
||||
}
|
||||
|
||||
BilinearForm a(&fespace);
|
||||
if (!vec)
|
||||
{
|
||||
a.AddDomainIntegrator(new MassIntegrator(oneCoef));
|
||||
}
|
||||
else
|
||||
{
|
||||
a.AddDomainIntegrator(new VectorFEMassIntegrator(oneCoef));
|
||||
}
|
||||
a.Assemble();
|
||||
|
||||
LinearForm Au(&fespace);
|
||||
LinearForm ATu(&fespace);
|
||||
LinearForm aAu(&fespace);
|
||||
LinearForm aATu(&fespace);
|
||||
LinearForm b(&fespace);
|
||||
|
||||
a.Mult(u, Au);
|
||||
a.MultTranspose(u, ATu);
|
||||
|
||||
aAu = Au;
|
||||
aATu = ATu;
|
||||
|
||||
a.AddMult(u, aAu, alpha);
|
||||
a.AddMultTranspose(u, aATu, alpha);
|
||||
|
||||
// Modify the Bilinear Form
|
||||
OperatorPtr A;
|
||||
a.FormSystemMatrix(ess_tdof_list, A);
|
||||
|
||||
a.FullMult(u, b);
|
||||
b -= Au;
|
||||
|
||||
REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
a.FullMultTranspose(u, b);
|
||||
b -= ATu;
|
||||
|
||||
REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
b = Au;
|
||||
a.FullAddMult(u, b, alpha);
|
||||
b -= aAu;
|
||||
|
||||
REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
b = ATu;
|
||||
a.FullAddMultTranspose(u, b, alpha);
|
||||
b -= aATu;
|
||||
|
||||
REQUIRE(b.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
delete fec;
|
||||
}
|
||||
}
|
||||
|
||||
delete mesh;
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("MixedBilinearform Full Ops",
|
||||
"[MixedBilinearForm]")
|
||||
{
|
||||
int order = 2;
|
||||
double alpha = M_E;
|
||||
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
mesh->UniformRefinement();
|
||||
|
||||
Vector oneVec(dim); oneVec = 1.0;
|
||||
|
||||
ConstantCoefficient oneCoef(1.0);
|
||||
VectorConstantCoefficient oneVecCoef(oneVec);
|
||||
|
||||
for (int ft = (int)FEType::H1_FEC; ft <= (int)FEType::RT_FEC; ft++)
|
||||
{
|
||||
bool vec = (ft == (int)FEType::ND_FEC || ft == (int)FEType::RT_FEC);
|
||||
|
||||
if (dim == 1 && vec) { continue; }
|
||||
if (vec && (mt == (int)MeshType::WEDGE2 ||
|
||||
mt == (int)MeshType::WEDGE4 ||
|
||||
mt == (int)MeshType::MIXED3D6 ||
|
||||
mt == (int)MeshType::MIXED3D8))
|
||||
{ continue; }
|
||||
|
||||
SECTION("Integral of field " + std::to_string(ft) +
|
||||
" on mesh type " + std::to_string(mt) )
|
||||
{
|
||||
|
||||
FiniteElementCollection *fec_dom = NULL;
|
||||
FiniteElementCollection *fec_ran = NULL;
|
||||
switch ((FEType)ft)
|
||||
{
|
||||
case FEType::H1_FEC:
|
||||
fec_dom = new H1_FECollection(order, dim);
|
||||
fec_ran = new H1_FECollection(order-1, dim);
|
||||
break;
|
||||
case FEType::ND_FEC:
|
||||
fec_dom = new ND_FECollection(order, dim);
|
||||
fec_ran = new RT_FECollection(order-1, dim);
|
||||
break;
|
||||
case FEType::RT_FEC:
|
||||
fec_dom = new RT_FECollection(order-1, dim);
|
||||
fec_ran = new ND_FECollection(order, dim);
|
||||
break;
|
||||
default:
|
||||
MFEM_ABORT("Invalid vector FE type");
|
||||
}
|
||||
FiniteElementSpace fespace_dom(mesh, fec_dom);
|
||||
FiniteElementSpace fespace_ran(mesh, fec_ran);
|
||||
|
||||
Array<int> ess_tdof_list_dom;
|
||||
Array<int> ess_tdof_list_ran;
|
||||
if (mesh->bdr_attributes.Size())
|
||||
{
|
||||
Array<int> ess_bdr(mesh->bdr_attributes.Max());
|
||||
ess_bdr = 1;
|
||||
fespace_dom.GetEssentialTrueDofs(ess_bdr, ess_tdof_list_dom);
|
||||
fespace_ran.GetEssentialTrueDofs(ess_bdr, ess_tdof_list_ran);
|
||||
}
|
||||
|
||||
GridFunction u_dom(&fespace_dom);
|
||||
GridFunction u_ran(&fespace_ran);
|
||||
if (!vec)
|
||||
{
|
||||
u_dom.ProjectCoefficient(oneCoef);
|
||||
u_ran.ProjectCoefficient(oneCoef);
|
||||
}
|
||||
else
|
||||
{
|
||||
u_dom.ProjectCoefficient(oneVecCoef);
|
||||
u_ran.ProjectCoefficient(oneVecCoef);
|
||||
}
|
||||
|
||||
MixedBilinearForm a(&fespace_dom, &fespace_ran);
|
||||
if (!vec)
|
||||
{
|
||||
a.AddDomainIntegrator(new MassIntegrator(oneCoef));
|
||||
}
|
||||
else
|
||||
{
|
||||
a.AddDomainIntegrator(new VectorFEMassIntegrator(oneCoef));
|
||||
}
|
||||
a.Assemble();
|
||||
|
||||
LinearForm Au(&fespace_ran);
|
||||
LinearForm ATu(&fespace_dom);
|
||||
LinearForm aAu(&fespace_ran);
|
||||
LinearForm aATu(&fespace_dom);
|
||||
LinearForm b_ran(&fespace_ran);
|
||||
LinearForm b_dom(&fespace_dom);
|
||||
|
||||
a.Mult(u_dom, Au);
|
||||
a.MultTranspose(u_ran, ATu);
|
||||
|
||||
aAu = Au;
|
||||
aATu = ATu;
|
||||
|
||||
a.AddMult(u_dom, aAu, alpha);
|
||||
a.AddMultTranspose(u_ran, aATu, alpha);
|
||||
|
||||
// Modify the Bilinear Form
|
||||
OperatorPtr A;
|
||||
a.FormRectangularSystemMatrix(ess_tdof_list_dom,
|
||||
ess_tdof_list_ran, A);
|
||||
|
||||
a.FullMult(u_dom, b_ran);
|
||||
b_ran -= Au;
|
||||
|
||||
REQUIRE(b_ran.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
a.FullMultTranspose(u_ran, b_dom);
|
||||
b_dom -= ATu;
|
||||
|
||||
REQUIRE(b_dom.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
b_ran = Au;
|
||||
a.FullAddMult(u_dom, b_ran, alpha);
|
||||
b_ran -= aAu;
|
||||
|
||||
REQUIRE(b_ran.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
b_dom = ATu;
|
||||
a.FullAddMultTranspose(u_ran, b_dom, alpha);
|
||||
b_dom -= aATu;
|
||||
|
||||
REQUIRE(b_dom.Norml2() == MFEM_Approx( 0.0));
|
||||
|
||||
delete fec_dom;
|
||||
delete fec_ran;
|
||||
}
|
||||
}
|
||||
|
||||
delete mesh;
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace bilinearform
|
||||
|
||||
@@ -11,6 +11,9 @@
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "unit_tests.hpp"
|
||||
#include "common_get_mesh.hpp"
|
||||
|
||||
using namespace mfem_test_fem;
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
@@ -22,29 +25,6 @@ static double a_ = M_PI;
|
||||
static double b_ = M_PI / sqrt(2.0);
|
||||
static double c_ = M_PI / 2.0;
|
||||
|
||||
enum MeshType
|
||||
{
|
||||
SEGMENT = 0,
|
||||
QUADRILATERAL = 1,
|
||||
TRIANGLE2A = 2,
|
||||
TRIANGLE2B = 3,
|
||||
TRIANGLE2C = 4,
|
||||
TRIANGLE4 = 5,
|
||||
MIXED2D = 6,
|
||||
HEXAHEDRON = 7,
|
||||
HEXAHEDRON2A = 8,
|
||||
HEXAHEDRON2B = 9,
|
||||
HEXAHEDRON2C = 10,
|
||||
HEXAHEDRON2D = 11,
|
||||
WEDGE2 = 12,
|
||||
TETRAHEDRA = 13,
|
||||
WEDGE4 = 14,
|
||||
MIXED3D6 = 15,
|
||||
MIXED3D8 = 16
|
||||
};
|
||||
|
||||
Mesh * GetMesh(MeshType type);
|
||||
|
||||
enum BasisType
|
||||
{
|
||||
H1 = 0, ND = 1, RT = 2, L2 = 3
|
||||
@@ -61,7 +41,7 @@ TEST_CASE("Build Dof To Arrays",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
if (dim < 3 ||
|
||||
mt == MeshType::HEXAHEDRON ||
|
||||
@@ -156,7 +136,7 @@ TEST_CASE("Build Dof To Arrays (Parallel)",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
if (dim < 3 ||
|
||||
mt == MeshType::HEXAHEDRON ||
|
||||
@@ -238,239 +218,6 @@ TEST_CASE("Build Dof To Arrays (Parallel)",
|
||||
}
|
||||
#endif // MFEM_USE_MPI
|
||||
|
||||
Mesh * GetMesh(MeshType type)
|
||||
{
|
||||
Mesh * mesh = NULL;
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case SEGMENT:
|
||||
mesh = new Mesh(1, 2, 1);
|
||||
mesh->AddVertex(0.0);
|
||||
mesh->AddVertex(a_);
|
||||
|
||||
mesh->AddSegment(0, 1);
|
||||
|
||||
mesh->AddBdrPoint(0);
|
||||
mesh->AddBdrPoint(1);
|
||||
break;
|
||||
case QUADRILATERAL:
|
||||
mesh = new Mesh(2, 4, 1);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddQuad(0, 1, 2, 3);
|
||||
break;
|
||||
case TRIANGLE2A:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddTriangle(0, 1, 2);
|
||||
mesh->AddTriangle(2, 3, 0);
|
||||
break;
|
||||
case TRIANGLE2B:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddTriangle(1, 2, 0);
|
||||
mesh->AddTriangle(3, 0, 2);
|
||||
break;
|
||||
case TRIANGLE2C:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddTriangle(2, 0, 1);
|
||||
mesh->AddTriangle(0, 2, 3);
|
||||
break;
|
||||
case TRIANGLE4:
|
||||
mesh = new Mesh(2, 5, 4);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
mesh->AddVertex(0.5 * a_, 0.5 * b_);
|
||||
|
||||
mesh->AddTriangle(0, 1, 4);
|
||||
mesh->AddTriangle(1, 2, 4);
|
||||
mesh->AddTriangle(2, 3, 4);
|
||||
mesh->AddTriangle(3, 0, 4);
|
||||
break;
|
||||
case MIXED2D:
|
||||
mesh = new Mesh(2, 6, 4);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
mesh->AddVertex(0.5 * b_, 0.5 * b_);
|
||||
mesh->AddVertex(a_ - 0.5 * b_, 0.5 * b_);
|
||||
|
||||
mesh->AddQuad(0, 1, 5, 4);
|
||||
mesh->AddTriangle(1, 2, 5);
|
||||
mesh->AddQuad(2, 3, 4, 5);
|
||||
mesh->AddTriangle(3, 0, 4);
|
||||
break;
|
||||
case HEXAHEDRON:
|
||||
mesh = new Mesh(3, 8, 1);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
|
||||
break;
|
||||
case HEXAHEDRON2A:
|
||||
case HEXAHEDRON2B:
|
||||
case HEXAHEDRON2C:
|
||||
case HEXAHEDRON2D:
|
||||
mesh = new Mesh(3, 12, 2);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(0.5 * a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.5 * a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(0.5 * a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.5 * a_, b_, c_);
|
||||
mesh->AddVertex(0.0,b_, c_);
|
||||
|
||||
mesh->AddHex(0, 5, 11, 6, 1, 4, 10, 7);
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case HEXAHEDRON2A: // Face Orientation 1
|
||||
mesh->AddHex(4, 10, 7, 1, 3, 9, 8, 2);
|
||||
break;
|
||||
case HEXAHEDRON2B: // Face Orientation 3
|
||||
mesh->AddHex(10, 7, 1, 4, 9, 8, 2, 3);
|
||||
break;
|
||||
case HEXAHEDRON2C: // Face Orientation 5
|
||||
mesh->AddHex(7, 1, 4, 10, 8, 2, 3, 9);
|
||||
break;
|
||||
case HEXAHEDRON2D: // Face Orientation 7
|
||||
mesh->AddHex(1, 4, 10, 7, 2, 3, 9, 8);
|
||||
break;
|
||||
default:
|
||||
// Cannot happen
|
||||
break;
|
||||
}
|
||||
break;
|
||||
case WEDGE2:
|
||||
mesh = new Mesh(3, 8, 2);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddWedge(0, 1, 2, 4, 5, 6);
|
||||
mesh->AddWedge(0, 2, 3, 4, 6, 7);
|
||||
break;
|
||||
case TETRAHEDRA:
|
||||
mesh = new Mesh(3, 8, 5);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddTet(0, 2, 7, 5);
|
||||
mesh->AddTet(6, 7, 2, 5);
|
||||
mesh->AddTet(4, 7, 5, 0);
|
||||
mesh->AddTet(1, 0, 5, 2);
|
||||
mesh->AddTet(3, 7, 0, 2);
|
||||
break;
|
||||
case WEDGE4:
|
||||
mesh = new Mesh(3, 10, 4);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.5 * a_, 0.5 * b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
mesh->AddVertex(0.5 * a_, 0.5 * b_, c_);
|
||||
|
||||
mesh->AddWedge(0, 1, 4, 5, 6, 9);
|
||||
mesh->AddWedge(1, 2, 4, 6, 7, 9);
|
||||
mesh->AddWedge(2, 3, 4, 7, 8, 9);
|
||||
mesh->AddWedge(3, 0, 4, 8, 5, 9);
|
||||
break;
|
||||
case MIXED3D6:
|
||||
mesh = new Mesh(3, 12, 6);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.5 * c_, 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(a_ - 0.5 * c_, 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(a_ - 0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
|
||||
mesh->AddWedge(0, 4, 8, 1, 5, 9);
|
||||
mesh->AddWedge(1, 5, 9, 2, 6, 10);
|
||||
mesh->AddWedge(2, 6, 10, 3, 7, 11);
|
||||
mesh->AddWedge(3, 7, 11, 0, 4, 8);
|
||||
mesh->AddHex(4, 5, 6, 7, 8, 9, 10, 11);
|
||||
break;
|
||||
case MIXED3D8:
|
||||
mesh = new Mesh(3, 10, 8);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.25 * a_, 0.5 * b_, 0.5 * c_);
|
||||
mesh->AddVertex(0.75 * a_, 0.5 * b_, 0.5 * c_);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddWedge(0, 3, 4, 1, 2, 5);
|
||||
mesh->AddWedge(3, 9, 4, 2, 8, 5);
|
||||
mesh->AddWedge(9, 6, 4, 8, 7, 5);
|
||||
mesh->AddWedge(6, 0, 4, 7, 1, 5);
|
||||
mesh->AddTet(0, 3, 9, 4);
|
||||
mesh->AddTet(0, 9, 6, 4);
|
||||
mesh->AddTet(1, 7, 2, 5);
|
||||
mesh->AddTet(8, 2, 7, 5);
|
||||
break;
|
||||
}
|
||||
mesh->FinalizeTopology();
|
||||
|
||||
return mesh;
|
||||
}
|
||||
|
||||
} // namespace build_dof_to_arrays
|
||||
|
||||
} // namespace mfem
|
||||
|
||||
@@ -11,8 +11,10 @@
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "unit_tests.hpp"
|
||||
#include "common_get_mesh.hpp"
|
||||
|
||||
using namespace mfem;
|
||||
using namespace mfem_test_fem;
|
||||
|
||||
namespace domain_int
|
||||
{
|
||||
@@ -46,29 +48,6 @@ enum FEType
|
||||
L2I_FEC,
|
||||
};
|
||||
|
||||
enum MeshType
|
||||
{
|
||||
SEGMENT = 0,
|
||||
QUADRILATERAL = 1,
|
||||
TRIANGLE2A = 2,
|
||||
TRIANGLE2B = 3,
|
||||
TRIANGLE2C = 4,
|
||||
TRIANGLE4 = 5,
|
||||
MIXED2D = 6,
|
||||
HEXAHEDRON = 7,
|
||||
HEXAHEDRON2A = 8,
|
||||
HEXAHEDRON2B = 9,
|
||||
HEXAHEDRON2C = 10,
|
||||
HEXAHEDRON2D = 11,
|
||||
WEDGE2 = 12,
|
||||
TETRAHEDRA = 13,
|
||||
WEDGE4 = 14,
|
||||
MIXED3D6 = 15,
|
||||
MIXED3D8 = 16
|
||||
};
|
||||
|
||||
Mesh * GetMesh(MeshType type);
|
||||
|
||||
TEST_CASE("Domain Integration (Scalar Field)",
|
||||
"[H1_FECollection]"
|
||||
"[L2_FECollection]"
|
||||
@@ -80,7 +59,7 @@ TEST_CASE("Domain Integration (Scalar Field)",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
mesh->UniformRefinement();
|
||||
|
||||
@@ -155,7 +134,7 @@ TEST_CASE("Domain Integration (Vector Field)",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
int sdim = mesh->SpaceDimension();
|
||||
mesh->UniformRefinement();
|
||||
@@ -258,7 +237,7 @@ TEST_CASE("Domain Integration in Parallel (Scalar Field)",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
while (mesh->GetNE() < num_procs)
|
||||
{
|
||||
@@ -340,7 +319,7 @@ TEST_CASE("Domain Integration in Parallel (Vector Field)",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
int sdim = mesh->SpaceDimension();
|
||||
while (mesh->GetNE() < num_procs)
|
||||
@@ -427,237 +406,4 @@ TEST_CASE("Domain Integration in Parallel (Vector Field)",
|
||||
|
||||
#endif // MFEM_USE_MPI
|
||||
|
||||
Mesh * GetMesh(MeshType type)
|
||||
{
|
||||
Mesh * mesh = NULL;
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case SEGMENT:
|
||||
mesh = new Mesh(1, 2, 1);
|
||||
mesh->AddVertex(0.0);
|
||||
mesh->AddVertex(a_);
|
||||
|
||||
mesh->AddSegment(0, 1);
|
||||
|
||||
mesh->AddBdrPoint(0);
|
||||
mesh->AddBdrPoint(1);
|
||||
break;
|
||||
case QUADRILATERAL:
|
||||
mesh = new Mesh(2, 4, 1);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddQuad(0, 1, 2, 3);
|
||||
break;
|
||||
case TRIANGLE2A:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddTriangle(0, 1, 2);
|
||||
mesh->AddTriangle(2, 3, 0);
|
||||
break;
|
||||
case TRIANGLE2B:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddTriangle(1, 2, 0);
|
||||
mesh->AddTriangle(3, 0, 2);
|
||||
break;
|
||||
case TRIANGLE2C:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
|
||||
mesh->AddTriangle(2, 0, 1);
|
||||
mesh->AddTriangle(0, 2, 3);
|
||||
break;
|
||||
case TRIANGLE4:
|
||||
mesh = new Mesh(2, 5, 4);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
mesh->AddVertex(0.5 * a_, 0.5 * b_);
|
||||
|
||||
mesh->AddTriangle(0, 1, 4);
|
||||
mesh->AddTriangle(1, 2, 4);
|
||||
mesh->AddTriangle(2, 3, 4);
|
||||
mesh->AddTriangle(3, 0, 4);
|
||||
break;
|
||||
case MIXED2D:
|
||||
mesh = new Mesh(2, 6, 4);
|
||||
mesh->AddVertex(0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0);
|
||||
mesh->AddVertex(a_, b_);
|
||||
mesh->AddVertex(0.0, b_);
|
||||
mesh->AddVertex(0.5 * b_, 0.5 * b_);
|
||||
mesh->AddVertex(a_ - 0.5 * b_, 0.5 * b_);
|
||||
|
||||
mesh->AddQuad(0, 1, 5, 4);
|
||||
mesh->AddTriangle(1, 2, 5);
|
||||
mesh->AddQuad(2, 3, 4, 5);
|
||||
mesh->AddTriangle(3, 0, 4);
|
||||
break;
|
||||
case HEXAHEDRON:
|
||||
mesh = new Mesh(3, 8, 1);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
|
||||
break;
|
||||
case HEXAHEDRON2A:
|
||||
case HEXAHEDRON2B:
|
||||
case HEXAHEDRON2C:
|
||||
case HEXAHEDRON2D:
|
||||
mesh = new Mesh(3, 12, 2);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(0.5 * a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.5 * a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(0.5 * a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.5 * a_, b_, c_);
|
||||
mesh->AddVertex(0.0,b_, c_);
|
||||
|
||||
mesh->AddHex(0, 5, 11, 6, 1, 4, 10, 7);
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case HEXAHEDRON2A: // Face Orientation 1
|
||||
mesh->AddHex(4, 10, 7, 1, 3, 9, 8, 2);
|
||||
break;
|
||||
case HEXAHEDRON2B: // Face Orientation 3
|
||||
mesh->AddHex(10, 7, 1, 4, 9, 8, 2, 3);
|
||||
break;
|
||||
case HEXAHEDRON2C: // Face Orientation 5
|
||||
mesh->AddHex(7, 1, 4, 10, 8, 2, 3, 9);
|
||||
break;
|
||||
case HEXAHEDRON2D: // Face Orientation 7
|
||||
mesh->AddHex(1, 4, 10, 7, 2, 3, 9, 8);
|
||||
break;
|
||||
default:
|
||||
// Cannot happen
|
||||
break;
|
||||
}
|
||||
break;
|
||||
case WEDGE2:
|
||||
mesh = new Mesh(3, 8, 2);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddWedge(0, 1, 2, 4, 5, 6);
|
||||
mesh->AddWedge(0, 2, 3, 4, 6, 7);
|
||||
break;
|
||||
case TETRAHEDRA:
|
||||
mesh = new Mesh(3, 8, 5);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddTet(0, 2, 7, 5);
|
||||
mesh->AddTet(6, 7, 2, 5);
|
||||
mesh->AddTet(4, 7, 5, 0);
|
||||
mesh->AddTet(1, 0, 5, 2);
|
||||
mesh->AddTet(3, 7, 0, 2);
|
||||
break;
|
||||
case WEDGE4:
|
||||
mesh = new Mesh(3, 10, 4);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.5 * a_, 0.5 * b_, 0.0);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
mesh->AddVertex(0.5 * a_, 0.5 * b_, c_);
|
||||
|
||||
mesh->AddWedge(0, 1, 4, 5, 6, 9);
|
||||
mesh->AddWedge(1, 2, 4, 6, 7, 9);
|
||||
mesh->AddWedge(2, 3, 4, 7, 8, 9);
|
||||
mesh->AddWedge(3, 0, 4, 8, 5, 9);
|
||||
break;
|
||||
case MIXED3D6:
|
||||
mesh = new Mesh(3, 12, 6);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.5 * c_, 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(a_ - 0.5 * c_, 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(a_ - 0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
|
||||
mesh->AddWedge(0, 4, 8, 1, 5, 9);
|
||||
mesh->AddWedge(1, 5, 9, 2, 6, 10);
|
||||
mesh->AddWedge(2, 6, 10, 3, 7, 11);
|
||||
mesh->AddWedge(3, 7, 11, 0, 4, 8);
|
||||
mesh->AddHex(4, 5, 6, 7, 8, 9, 10, 11);
|
||||
break;
|
||||
case MIXED3D8:
|
||||
mesh = new Mesh(3, 10, 8);
|
||||
mesh->AddVertex(0.0, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, 0.0, 0.0);
|
||||
mesh->AddVertex(a_, b_, 0.0);
|
||||
mesh->AddVertex(0.0, b_, 0.0);
|
||||
mesh->AddVertex(0.25 * a_, 0.5 * b_, 0.5 * c_);
|
||||
mesh->AddVertex(0.75 * a_, 0.5 * b_, 0.5 * c_);
|
||||
mesh->AddVertex(0.0, 0.0, c_);
|
||||
mesh->AddVertex(a_, 0.0, c_);
|
||||
mesh->AddVertex(a_, b_, c_);
|
||||
mesh->AddVertex(0.0, b_, c_);
|
||||
|
||||
mesh->AddWedge(0, 3, 4, 1, 2, 5);
|
||||
mesh->AddWedge(3, 9, 4, 2, 8, 5);
|
||||
mesh->AddWedge(9, 6, 4, 8, 7, 5);
|
||||
mesh->AddWedge(6, 0, 4, 7, 1, 5);
|
||||
mesh->AddTet(0, 3, 9, 4);
|
||||
mesh->AddTet(0, 9, 6, 4);
|
||||
mesh->AddTet(1, 7, 2, 5);
|
||||
mesh->AddTet(8, 2, 7, 5);
|
||||
break;
|
||||
}
|
||||
mesh->FinalizeTopology();
|
||||
|
||||
return mesh;
|
||||
}
|
||||
|
||||
} // namespace domain_int
|
||||
|
||||
@@ -11,38 +11,21 @@
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include "unit_tests.hpp"
|
||||
#include "common_get_mesh.hpp"
|
||||
|
||||
using namespace mfem;
|
||||
using namespace mfem_test_fem;
|
||||
|
||||
namespace eigs
|
||||
{
|
||||
|
||||
#if defined MFEM_USE_LAPACK || defined MFEM_USE_MPI
|
||||
|
||||
static double a_ = M_PI;
|
||||
static double b_ = M_PI / sqrt(2.0);
|
||||
static double c_ = M_PI / 2.0;
|
||||
|
||||
enum MeshType
|
||||
{
|
||||
SEGMENT = 0,
|
||||
QUADRILATERAL = 1,
|
||||
TRIANGLE2A = 2,
|
||||
TRIANGLE2B = 3,
|
||||
TRIANGLE2C = 4,
|
||||
TRIANGLE4 = 5,
|
||||
MIXED2D = 6,
|
||||
HEXAHEDRON = 7,
|
||||
HEXAHEDRON2A = 8,
|
||||
HEXAHEDRON2B = 9,
|
||||
HEXAHEDRON2C = 10,
|
||||
HEXAHEDRON2D = 11,
|
||||
WEDGE2 = 12,
|
||||
TETRAHEDRA = 13,
|
||||
WEDGE4 = 14,
|
||||
MIXED3D6 = 15,
|
||||
MIXED3D8 = 16
|
||||
};
|
||||
|
||||
Mesh * GetMesh(MeshType type);
|
||||
#endif
|
||||
|
||||
int eigs[21] =
|
||||
{
|
||||
@@ -63,7 +46,7 @@ TEST_CASE("Laplacian Eigenvalues",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
if (dim < 3 ||
|
||||
mt == MeshType::HEXAHEDRON ||
|
||||
@@ -178,7 +161,7 @@ TEST_CASE("Laplacian Eigenvalues in Parallel",
|
||||
for (int mt = (int)MeshType::SEGMENT;
|
||||
mt <= (int)MeshType::MIXED3D8; mt++)
|
||||
{
|
||||
Mesh *mesh = GetMesh((MeshType)mt);
|
||||
Mesh *mesh = GetMesh((MeshType)mt, a_, b_, c_);
|
||||
int dim = mesh->Dimension();
|
||||
if (dim < 3 ||
|
||||
mt == MeshType::HEXAHEDRON ||
|
||||
@@ -273,401 +256,4 @@ TEST_CASE("Laplacian Eigenvalues in Parallel",
|
||||
|
||||
#endif // MFEM_USE_MPI
|
||||
|
||||
Mesh * GetMesh(MeshType type)
|
||||
{
|
||||
Mesh * mesh = NULL;
|
||||
double c[3];
|
||||
int v[8];
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case SEGMENT:
|
||||
mesh = new Mesh(1, 2, 1);
|
||||
c[0] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_;
|
||||
mesh->AddVertex(c);
|
||||
v[0] = 0; v[1] = 1;
|
||||
mesh->AddSegment(v);
|
||||
{
|
||||
Element * el = mesh->NewElement(Geometry::POINT);
|
||||
el->SetAttribute(1);
|
||||
el->SetVertices(&v[0]);
|
||||
mesh->AddBdrElement(el);
|
||||
}
|
||||
{
|
||||
Element * el = mesh->NewElement(Geometry::POINT);
|
||||
el->SetAttribute(2);
|
||||
el->SetVertices(&v[1]);
|
||||
mesh->AddBdrElement(el);
|
||||
}
|
||||
break;
|
||||
case QUADRILATERAL:
|
||||
mesh = new Mesh(2, 4, 1);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
|
||||
mesh->AddQuad(v);
|
||||
break;
|
||||
case TRIANGLE2A:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 0;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case TRIANGLE2B:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 1; v[1] = 2; v[2] = 0;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 2;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case TRIANGLE2C:
|
||||
mesh = new Mesh(2, 4, 2);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 2; v[1] = 0; v[2] = 1;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 0; v[1] = 2; v[2] = 3;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case TRIANGLE4:
|
||||
mesh = new Mesh(2, 5, 4);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = 0.5 * b_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 1; v[1] = 2; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case MIXED2D:
|
||||
mesh = new Mesh(2, 6, 4);
|
||||
c[0] = 0.0; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * b_; c[1] = 0.5 * b_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_ - 0.5 * b_; c[1] = 0.5 * b_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 5; v[3] = 4;
|
||||
mesh->AddQuad(v);
|
||||
v[0] = 1; v[1] = 2; v[2] = 5;
|
||||
mesh->AddTri(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 5;
|
||||
mesh->AddQuad(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 4;
|
||||
mesh->AddTri(v);
|
||||
break;
|
||||
case HEXAHEDRON:
|
||||
mesh = new Mesh(3, 8, 1);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
|
||||
v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2A:
|
||||
case HEXAHEDRON2B:
|
||||
case HEXAHEDRON2C:
|
||||
case HEXAHEDRON2D:
|
||||
mesh = new Mesh(3, 12, 2);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 5; v[2] = 11; v[3] = 6;
|
||||
v[4] = 1; v[5] = 4; v[6] = 10; v[7] = 7;
|
||||
mesh->AddHex(v);
|
||||
|
||||
switch (type)
|
||||
{
|
||||
case HEXAHEDRON2A: // Face Orientation 1
|
||||
v[0] = 4; v[1] = 10; v[2] = 7; v[3] = 1;
|
||||
v[4] = 3; v[5] = 9; v[6] = 8; v[7] = 2;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2B: // Face Orientation 3
|
||||
v[0] = 10; v[1] = 7; v[2] = 1; v[3] = 4;
|
||||
v[4] = 9; v[5] = 8; v[6] = 2; v[7] = 3;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2C: // Face Orientation 5
|
||||
v[0] = 7; v[1] = 1; v[2] = 4; v[3] = 10;
|
||||
v[4] = 8; v[5] = 2; v[6] = 3; v[7] = 9;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case HEXAHEDRON2D: // Face Orientation 7
|
||||
v[0] = 1; v[1] = 4; v[2] = 10; v[3] = 7;
|
||||
v[4] = 2; v[5] = 3; v[6] = 9; v[7] = 8;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
default:
|
||||
// Cannot happen
|
||||
break;
|
||||
}
|
||||
break;
|
||||
case WEDGE2:
|
||||
mesh = new Mesh(3, 8, 2);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 4; v[4] = 5; v[5] = 6;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 0; v[1] = 2; v[2] = 3; v[3] = 4; v[4] = 6; v[5] = 7;
|
||||
mesh->AddWedge(v);
|
||||
break;
|
||||
case TETRAHEDRA:
|
||||
mesh = new Mesh(3, 8, 5);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 2; v[2] = 7; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 6; v[1] = 7; v[2] = 2; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 4; v[1] = 7; v[2] = 5; v[3] = 0;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 1; v[1] = 0; v[2] = 5; v[3] = 2;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 3; v[1] = 7; v[2] = 0; v[3] = 2;
|
||||
mesh->AddTet(v);
|
||||
break;
|
||||
case WEDGE4:
|
||||
mesh = new Mesh(3, 10, 4);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = 0.5 * b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * a_; c[1] = 0.5 * b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 4; v[3] = 5; v[4] = 6; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 1; v[1] = 2; v[2] = 4; v[3] = 6; v[4] = 7; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 2; v[1] = 3; v[2] = 4; v[3] = 7; v[4] = 8; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 3; v[1] = 0; v[2] = 4; v[3] = 8; v[4] = 5; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
break;
|
||||
case MIXED3D6:
|
||||
mesh = new Mesh(3, 12, 6);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * c_; c[1] = 0.5 * c_; c[2] = 0.5 * c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_ - 0.5 * c_; c[1] = 0.5 * c_; c[2] = 0.5 * c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_ - 0.5 * c_; c[1] = b_ - 0.5 * c_; c[2] = 0.5 * c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.5 * c_; c[1] = b_ - 0.5 * c_; c[2] = 0.5 * c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3;
|
||||
v[4] = 4; v[5] = 5; v[6] = 6; v[7] = 7;
|
||||
mesh->AddHex(v);
|
||||
v[0] = 0; v[1] = 4; v[2] = 8; v[3] = 1; v[4] = 5; v[5] = 9;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 1; v[1] = 5; v[2] = 9; v[3] = 2; v[4] = 6; v[5] = 10;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 2; v[1] = 6; v[2] = 10; v[3] = 3; v[4] = 7; v[5] = 11;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 3; v[1] = 7; v[2] = 11; v[3] = 0; v[4] = 4; v[5] = 8;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 4; v[1] = 5; v[2] = 6; v[3] = 7;
|
||||
v[4] = 8; v[5] = 9; v[6] = 10; v[7] = 11;
|
||||
mesh->AddHex(v);
|
||||
break;
|
||||
case MIXED3D8:
|
||||
mesh = new Mesh(3, 10, 8);
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = 0.0;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
c[0] = 0.25 * a_; c[1] = 0.5 * b_; c[2] = 0.5 * c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.75 * a_; c[1] = 0.5 * b_; c[2] = 0.5 * c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
c[0] = 0.0; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = 0.0; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = a_; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
c[0] = 0.0; c[1] = b_; c[2] = c_;
|
||||
mesh->AddVertex(c);
|
||||
|
||||
v[0] = 0; v[1] = 3; v[2] = 4; v[3] = 1; v[4] = 2; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 3; v[1] = 9; v[2] = 4; v[3] = 2; v[4] = 8; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 9; v[1] = 6; v[2] = 4; v[3] = 8; v[4] = 7; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 6; v[1] = 0; v[2] = 4; v[3] = 7; v[4] = 1; v[5] = 5;
|
||||
mesh->AddWedge(v);
|
||||
v[0] = 0; v[1] = 3; v[2] = 9; v[3] = 4;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 0; v[1] = 9; v[2] = 6; v[3] = 4;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 1; v[1] = 7; v[2] = 2; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
v[0] = 8; v[1] = 2; v[2] = 7; v[3] = 5;
|
||||
mesh->AddTet(v);
|
||||
break;
|
||||
}
|
||||
mesh->FinalizeTopology();
|
||||
|
||||
return mesh;
|
||||
}
|
||||
|
||||
} // namespace eigs
|
||||
|
||||
Reference in New Issue
Block a user