MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
toroid.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// Toroid Miniapp: Generate simple toroidal meshes
14// ------------------------------------------------
15//
16// This miniapp generates two types of toroidal meshes; one with triangular
17// cross sections and one with square cross sections. It works by defining a
18// stack of individual elements and bending them so that the bottom and top of
19// the stack can be joined to form a torus. The stack can also be twisted so
20// that the vertices of the bottom and top can be joined with any integer
21// offset.
22//
23// Compile with: make toroid
24//
25// Sample runs: toroid
26// toroid -nphi 6
27// toroid -ns 1
28// toroid -ns 0 -t0 -30
29// toroid -R 2 -r 1 -ns 3
30// toroid -R 2 -r 1 -ns -3
31// toroid -R 2 -r 1 -ns 3 -e 1
32// toroid -R 2 -r 1 -ns 3 -e 1 -rs 1
33// toroid -nphi 2 -ns 10 -e 1 -o 4
34
35#include "mfem.hpp"
36#include <fstream>
37#include <iostream>
38
39using namespace std;
40using namespace mfem;
41
42static Element::Type el_type_ = Element::WEDGE;
43static int order_ = 3;
44static int nphi_ = 8;
45static int ns_ = 0;
46static real_t R_ = 1.0;
47static real_t r_ = 0.2;
48static real_t theta0_ = 0.0;
49
50void pts(int iphi, int t, real_t x[]);
51void trans(const Vector &x, Vector &p);
52
53int main(int argc, char *argv[])
54{
55 int ser_ref_levels = 0;
56 int el_type = 0;
57 bool dg_mesh = false;
58 bool visualization = true;
59
60 OptionsParser args(argc, argv);
62 "Number of elements in phi-direction.");
64 "Number of shifts.");
66 "Order (polynomial degree) of the mesh elements.");
68 "Number of times to refine the mesh uniformly in serial.");
70 "Major radius of the torus.");
72 "Minor radius of the torus.");
74 "Starting angle of the cross section (in degrees).");
76 "Element type: 0 - Wedge, 1 - Hexahedron.");
77 args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
78 "Use discontinuous or continuous space for the mesh nodes.");
80 "--no-visualization",
81 "Enable or disable GLVis visualization.");
82 args.Parse();
83 if (!args.Good())
84 {
85 args.PrintUsage(cout);
86 return 1;
87 }
88 args.PrintOptions(cout);
89
90 // The output mesh could be hexahedra or prisms
91 el_type_ = (el_type == 0) ? Element::WEDGE : Element::HEXAHEDRON;
92 if (el_type_ != Element::WEDGE && el_type_ != Element::HEXAHEDRON)
93 {
94 cout << "Unsupported element type" << endl;
95 exit(1);
96 }
97
98 // Determine the number of nodes in the cross section
99 int nnode = (el_type_ == Element::WEDGE)? 3:4;
100 int nshift = (ns_ >= 0) ? 0 : (nnode * (1 - ns_ / nnode));
101
102 // Convert initial angle from degrees to radians
103 theta0_ *= M_PI / 180.0;
104
105 // Define an empty mesh
106 Mesh *mesh;
107 mesh = new Mesh(3, nnode * (nphi_+1), nphi_);
108
109 // Add vertices for a stack of elements
110 real_t c[3];
111 for (int i=0; i<=nphi_; i++)
112 {
113 c[0] = 0.0; c[1] = 0.0; c[2] = i;
115
116 c[0] = 1.0;
118
119 if (el_type_ == Element::HEXAHEDRON)
120 {
121 c[0] = 1.0; c[1] = 1.0;
123 }
124
125 c[0] = 0.0; c[1] = 1.0;
127 }
128
129 // Add Elements of the desired type
130 {
131 int v[8];
132 for (int i=0; i < nphi_; i++)
133 {
134 if (el_type_ == Element::WEDGE)
135 {
136 for (int j = 0; j < 6; j++) { v[j] = 3*i+j; }
138 }
139 else
140 {
141 for (int j = 0; j < 8; j++) { v[j] = 4*i+j; }
143 }
144 }
145 }
146 mesh->FinalizeTopology();
147
148 // Promote to high order mesh and transform into a torus shape
149 if (order_ > 1)
150 {
151 mesh->SetCurvature(order_, true, 3, Ordering::byVDIM);
152 }
153 mesh->Transform(trans);
154
155 // Stitch the ends of the stack together
156 {
157 Array<int> v2v(mesh->GetNV());
158 for (int i = 0; i < v2v.Size() - nnode; i++)
159 {
160 v2v[i] = i;
161 }
162 // identify vertices at the extremes of the stack of prisms
163 for (int i=0; i<nnode; i++)
164 {
165 v2v[v2v.Size() - nnode + i] = (nshift + ns_ + i) % nnode;
166 }
167 // renumber elements
168 for (int i = 0; i < mesh->GetNE(); i++)
169 {
170 Element *el = mesh->GetElement(i);
171 int *v = el->GetVertices();
172 int nv = el->GetNVertices();
173 for (int j = 0; j < nv; j++)
174 {
175 v[j] = v2v[v[j]];
176 }
177 }
178 // renumber boundary elements
179 for (int i = 0; i < mesh->GetNBE(); i++)
180 {
181 Element *el = mesh->GetBdrElement(i);
182 int *v = el->GetVertices();
183 int nv = el->GetNVertices();
184 for (int j = 0; j < nv; j++)
185 {
186 v[j] = v2v[v[j]];
187 }
188 }
189 mesh->RemoveUnusedVertices();
191 }
192 if (order_ > 1)
193 {
194 mesh->SetCurvature(order_, dg_mesh, 3, Ordering::byVDIM);
195 }
196
197 // Refine the mesh if desired
198 for (int lev = 0; lev < ser_ref_levels; lev++)
199 {
200 mesh->UniformRefinement();
201 }
202
203 // Output the resulting mesh to a file
204 {
205 ostringstream oss;
206 if (el_type_ == Element::WEDGE)
207 {
208 oss << "toroid-wedge";
209 }
210 else
211 {
212 oss << "toroid-hex";
213 }
214 oss << "-o" << order_ << "-s" << ns_;
215 if (ser_ref_levels > 0)
216 {
217 oss << "-r" << ser_ref_levels;
218 }
219 oss << ".mesh";
220 ofstream ofs(oss.str().c_str());
221 ofs.precision(8);
222 mesh->Print(ofs);
223 ofs.close();
224 }
225
226 // Output the resulting mesh to GLVis
227 if (visualization)
228 {
229 char vishost[] = "localhost";
230 int visport = 19916;
231 socketstream sol_sock(vishost, visport);
232 sol_sock.precision(8);
233 sol_sock << "mesh\n" << *mesh << flush;
234 }
235
236 // Clean up and exit
237 delete mesh;
238 return 0;
239}
240
241void trans(const Vector &x, Vector &p)
242{
243 int nnode = (el_type_ == Element::WEDGE)? 3:4;
244
245 real_t phi = 2.0 * M_PI * x[2] / nphi_;
246 real_t theta = theta0_ + phi * ns_ / nnode;
247
248 real_t u = (1.5 * (x[0] + x[1]) - 1.0) * r_;
249 real_t v = sqrt(0.75) * (x[0] - x[1]) * r_;
250
251 if (el_type_ == Element::WEDGE)
252 {
253 u = (1.5 * (x[0] + x[1]) - 1.0) * r_;
254 v = sqrt(0.75) * (x[0] - x[1]) * r_;
255 }
256 else
257 {
258 u = M_SQRT2 * (x[1] - 0.5) * r_;
259 v = M_SQRT2 * (x[0] - 0.5) * r_;
260 }
261
262 p[0] = ( R_ + u * cos(theta) + v * sin(theta)) * cos(phi);
263 p[1] = ( R_ + u * cos(theta) + v * sin(theta)) * sin(phi);
264 p[2] = v * cos(theta) - u * sin(theta);
265}
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
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
Adds a wedge to the mesh given by 6 vertices v1 through v6.
Definition mesh.cpp:1778
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition mesh.cpp:3135
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
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
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 AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Adds a hexahedron to the mesh given by 8 vertices v1 through v8.
Definition mesh.cpp:1806
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
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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 toroid.cpp:241
void pts(int iphi, int t, real_t x[])