MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
twist.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// -------------------------------------------------------
13// Twist Miniapp: Generate simple twisted periodic meshes
14// -------------------------------------------------------
15//
16// This miniapp generates simple periodic meshes to demonstrate MFEM's handling
17// of periodic domains. MFEM's strategy is to use a discontinuous vector field
18// to define the mesh coordinates on a topologically periodic mesh. It works by
19// defining a stack of individual elements and stitching together the top and
20// bottom of the mesh. The stack can also be twisted so that the vertices of the
21// bottom and top can be joined with any integer offset (for tetrahedral and
22// wedge meshes only even offsets are supported).
23//
24// Compile with: make twist
25//
26// Sample runs: twist
27// twist -no-pm
28// twist -nt -2 -no-pm
29// twist -nt 2 -e 4
30// twist -nt 2 -e 6
31// twist -nt 3 -e 8
32//
33#include "mfem.hpp"
34#include <fstream>
35#include <iostream>
36
37using namespace std;
38using namespace mfem;
39
40static Element::Type el_type_ = Element::WEDGE;
41static int order_ = 3;
42static int nz_ = 3;
43static int nt_ = 2;
44static real_t a_ = 1.0;
45static real_t b_ = 1.0;
46static real_t c_ = 3.0;
47
48void pts(int iphi, int t, real_t x[]);
49void trans(const Vector &x, Vector &p);
50
51int main(int argc, char *argv[])
52{
53 int ser_ref_levels = 0;
54 int el_type = 8;
55 bool per_mesh = true;
56 bool dg_mesh = false;
57 bool visualization = true;
58
59 OptionsParser args(argc, argv);
61 "Number of elements in z-direction.");
63 "Number of node positions to twist the top of the mesh.");
65 "Order (polynomial degree) of the mesh elements.");
67 "Number of times to refine the mesh uniformly in serial.");
69 "Width of the base in x-direction.");
71 "Width of the base in y-direction.");
73 "Height in z-direction.");
75 "Element type: 4 - Tetrahedron, 6 - Wedge, 8 - Hexahedron.");
77 "-no-pm", "--non-periodic-mesh",
78 "Enforce periodicity in z-direction "
79 "(requires discontinuous mesh).");
80 args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
81 "Use discontinuous or continuous space for the mesh nodes.");
83 "--no-visualization",
84 "Enable or disable GLVis visualization.");
85 args.Parse();
86 if (!args.Good())
87 {
88 args.PrintUsage(cout);
89 return 1;
90 }
91 args.PrintOptions(cout);
92
93 // The output mesh could be tetrahedra, hexahedra, or prisms
94 switch (el_type)
95 {
96 case 4:
97 el_type_ = Element::TETRAHEDRON;
98 break;
99 case 6:
100 el_type_ = Element::WEDGE;
101 break;
102 case 8:
103 el_type_ = Element::HEXAHEDRON;
104 break;
105 default:
106 cout << "Unsupported element type" << endl;
107 exit(1);
108 break;
109 }
110
111 // Define the mesh
112 Mesh mesh = Mesh::MakeCartesian3D(1, 1, nz_, el_type_, a_, b_, c_, false);
113
114 // Promote to high order mesh and transform into a twisted shape
115 if (order_ > 1 || dg_mesh || per_mesh)
116 {
117 mesh.SetCurvature(order_, dg_mesh || per_mesh, 3, Ordering::byVDIM);
118 }
119 if (nt_ != 0 )
120 {
121 mesh.Transform(trans);
122 }
123
124 while (per_mesh)
125 {
126 // Verify geometric compatibility
127 if (nt_ % 2 == 1 && std::abs(a_ - b_) > 1e-6 * a_)
128 {
129 cout << "Base is rectangular so number of shifts must be even "
130 << "for a periodic mesh!" << endl;
131 exit(1);
132 }
133
134 // Verify topological compatibility
135 if (nt_ % 2 == 1 && (el_type_ == Element::TETRAHEDRON ||
136 el_type_ == Element::WEDGE))
137 {
138 cout << "Diagonal cuts on the base and top must line up "
139 << "for a periodic mesh!" << endl;
140 exit(1);
141 }
142
143 int nnode = 4;
144 int noff = (nt_ >= 0) ? 0 : (nnode * (1 - nt_ / nnode));
145
146 Array<int> v2v(mesh.GetNV());
147 for (int i = 0; i < v2v.Size() - nnode; i++)
148 {
149 v2v[i] = i;
150 }
151 // identify vertices at the extremes of the stack
152 switch ((noff + nt_) % nnode)
153 {
154 case 0:
155 v2v[v2v.Size() - nnode + 0] = 0;
156 v2v[v2v.Size() - nnode + 1] = 1;
157 v2v[v2v.Size() - nnode + 2] = 2;
158 v2v[v2v.Size() - nnode + 3] = 3;
159 break;
160 case 1:
161 v2v[v2v.Size() - nnode + 0] = 2;
162 v2v[v2v.Size() - nnode + 1] = 0;
163 v2v[v2v.Size() - nnode + 2] = 3;
164 v2v[v2v.Size() - nnode + 3] = 1;
165 break;
166 case 2:
167 v2v[v2v.Size() - nnode + 0] = 3;
168 v2v[v2v.Size() - nnode + 1] = 2;
169 v2v[v2v.Size() - nnode + 2] = 1;
170 v2v[v2v.Size() - nnode + 3] = 0;
171 break;
172 case 3:
173 v2v[v2v.Size() - nnode + 0] = 1;
174 v2v[v2v.Size() - nnode + 1] = 3;
175 v2v[v2v.Size() - nnode + 2] = 0;
176 v2v[v2v.Size() - nnode + 3] = 2;
177 break;
178 }
179 // renumber elements
180 for (int i = 0; i < mesh.GetNE(); i++)
181 {
182 Element *el = mesh.GetElement(i);
183 int *v = el->GetVertices();
184 int nv = el->GetNVertices();
185 for (int j = 0; j < nv; j++)
186 {
187 v[j] = v2v[v[j]];
188 }
189 }
190 // renumber boundary elements
191 for (int i = 0; i < mesh.GetNBE(); i++)
192 {
193 Element *el = mesh.GetBdrElement(i);
194 int *v = el->GetVertices();
195 int nv = el->GetNVertices();
196 for (int j = 0; j < nv; j++)
197 {
198 v[j] = v2v[v[j]];
199 }
200 }
203
204 break;
205 }
206
207 // Refine the mesh if desired
208 for (int lev = 0; lev < ser_ref_levels; lev++)
209 {
210 mesh.UniformRefinement();
211 }
212
213 // Output the resulting mesh to a file
214 {
215 ostringstream oss;
216 if (el_type_ == Element::TETRAHEDRON)
217 {
218 oss << "twist-tet";
219 }
220 else if (el_type_ == Element::WEDGE)
221 {
222 oss << "twist-wedge";
223 }
224 else
225 {
226 oss << "twist-hex";
227 }
228 oss << "-o" << order_ << "-s" << nt_;
229 if (ser_ref_levels > 0)
230 {
231 oss << "-r" << ser_ref_levels;
232 }
233 if (per_mesh)
234 {
235 oss << "-p";
236 }
237 else if (dg_mesh)
238 {
239 oss << "-d";
240 }
241 else
242 {
243 oss << "-c";
244 }
245
246 oss << ".mesh";
247 ofstream ofs(oss.str().c_str());
248 ofs.precision(8);
249 mesh.Print(ofs);
250 ofs.close();
251 }
252
253 // Output the resulting mesh to GLVis
254 if (visualization)
255 {
256 char vishost[] = "localhost";
257 int visport = 19916;
258 socketstream sol_sock(vishost, visport);
259 sol_sock.precision(8);
260 sol_sock << "mesh\n" << mesh << flush;
261 }
262
263 // Clean up and exit
264 return 0;
265}
266
267void trans(const Vector &x, Vector &p)
268{
269 real_t z = x[2];
270 real_t phi = 0.5 * M_PI * nt_ * z / c_;
271 real_t cp = cos(phi);
272 real_t sp = sin(phi);
273
274 p[0] = 0.5 * a_ + (x[0] - 0.5 * a_) * cp - (x[1] - 0.5 * b_) * sp;
275 p[1] = 0.5 * b_ + (x[0] - 0.5 * a_) * sp + (x[1] - 0.5 * b_) * cp;
276 p[2] = z;
277}
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
Type
Constants for the classes derived from Element.
Definition element.hpp:41
virtual int GetNVertices() const =0
Mesh data type.
Definition mesh.hpp:56
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0, bool sfc_ordering=true)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
Definition mesh.cpp:4253
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition mesh.hpp:1223
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition mesh.cpp:6211
void Transform(void(*f)(const Vector &, Vector &))
Definition mesh.cpp:12821
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
void RemoveInternalBoundaries()
Definition mesh.cpp:12979
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition mesh.cpp:12872
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Vector data type.
Definition vector.hpp:80
int main()
const int visport
float real_t
Definition config.hpp:43
const char vishost[]
real_t p(const Vector &x, real_t t)
RefCoord t[3]
void trans(const Vector &x, Vector &p)
Definition twist.cpp:267
void pts(int iphi, int t, real_t x[])