MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
toroid.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 // ------------------------------------------------
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 
39 using namespace std;
40 using namespace mfem;
41 
42 static Element::Type el_type_ = Element::WEDGE;
43 static int order_ = 3;
44 static int nphi_ = 8;
45 static int ns_ = 0;
46 static double R_ = 1.0;
47 static double r_ = 0.2;
48 static double theta0_ = 0.0;
49 
50 void pts(int iphi, int t, double x[]);
51 void trans(const Vector &x, Vector &p);
52 
53 int 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);
61  args.AddOption(&nphi_, "-nphi", "--num-elements-phi",
62  "Number of elements in phi-direction.");
63  args.AddOption(&ns_, "-ns", "--num-shifts",
64  "Number of shifts.");
65  args.AddOption(&order_, "-o", "--mesh-order",
66  "Order (polynomial degree) of the mesh elements.");
67  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
68  "Number of times to refine the mesh uniformly in serial.");
69  args.AddOption(&R_, "-R", "--major-radius",
70  "Major radius of the torus.");
71  args.AddOption(&r_, "-r", "--minor-radius",
72  "Minor radius of the torus.");
73  args.AddOption(&theta0_, "-t0", "--initial-angle",
74  "Starting angle of the cross section (in degrees).");
75  args.AddOption(&el_type, "-e", "--element-type",
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.");
79  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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  double c[3];
111  for (int i=0; i<=nphi_; i++)
112  {
113  c[0] = 0.0; c[1] = 0.0; c[2] = i;
114  mesh->AddVertex(c);
115 
116  c[0] = 1.0;
117  mesh->AddVertex(c);
118 
119  if (el_type_ == Element::HEXAHEDRON)
120  {
121  c[0] = 1.0; c[1] = 1.0;
122  mesh->AddVertex(c);
123  }
124 
125  c[0] = 0.0; c[1] = 1.0;
126  mesh->AddVertex(c);
127  }
128 
129  // Add Elements of the desired type
130  int v[8];
131  for (int i=0; i < nphi_; i++)
132  {
133  if (el_type_ == Element::WEDGE)
134  {
135  for (int j = 0; j < 6; j++) { v[j] = 3*i+j; }
136  mesh->AddWedge(v);
137  }
138  else
139  {
140  for (int j = 0; j < 8; j++) { v[j] = 4*i+j; }
141  mesh->AddHex(v);
142  }
143  }
144  mesh->FinalizeTopology();
145 
146  // Promote to high order mesh and transform into a torus shape
147  if (order_ > 1)
148  {
149  mesh->SetCurvature(order_, true, 3, Ordering::byVDIM);
150  }
151  mesh->Transform(trans);
152 
153  // Stitch the ends of the stack together
154  {
155  Array<int> v2v(mesh->GetNV());
156  for (int i = 0; i < v2v.Size() - nnode; i++)
157  {
158  v2v[i] = i;
159  }
160  // identify vertices at the extremes of the stack of prisms
161  for (int i=0; i<nnode; i++)
162  {
163  v2v[v2v.Size() - nnode + i] = (nshift + ns_ + i) % nnode;
164  }
165  // renumber elements
166  for (int i = 0; i < mesh->GetNE(); i++)
167  {
168  Element *el = mesh->GetElement(i);
169  int *v = el->GetVertices();
170  int nv = el->GetNVertices();
171  for (int j = 0; j < nv; j++)
172  {
173  v[j] = v2v[v[j]];
174  }
175  }
176  // renumber boundary elements
177  for (int i = 0; i < mesh->GetNBE(); i++)
178  {
179  Element *el = mesh->GetBdrElement(i);
180  int *v = el->GetVertices();
181  int nv = el->GetNVertices();
182  for (int j = 0; j < nv; j++)
183  {
184  v[j] = v2v[v[j]];
185  }
186  }
187  mesh->RemoveUnusedVertices();
188  mesh->RemoveInternalBoundaries();
189  }
190  if (order_ > 1)
191  {
192  mesh->SetCurvature(order_, dg_mesh, 3, Ordering::byVDIM);
193  }
194 
195  // Refine the mesh if desired
196  for (int lev = 0; lev < ser_ref_levels; lev++)
197  {
198  mesh->UniformRefinement();
199  }
200 
201  // Output the resulting mesh to a file
202  {
203  ostringstream oss;
204  if (el_type_ == Element::WEDGE)
205  {
206  oss << "toroid-wedge";
207  }
208  else
209  {
210  oss << "toroid-hex";
211  }
212  oss << "-o" << order_ << "-s" << ns_;
213  if (ser_ref_levels > 0)
214  {
215  oss << "-r" << ser_ref_levels;
216  }
217  oss << ".mesh";
218  ofstream ofs(oss.str().c_str());
219  ofs.precision(8);
220  mesh->Print(ofs);
221  ofs.close();
222  }
223 
224  // Output the resulting mesh to GLVis
225  if (visualization)
226  {
227  char vishost[] = "localhost";
228  int visport = 19916;
229  socketstream sol_sock(vishost, visport);
230  sol_sock.precision(8);
231  sol_sock << "mesh\n" << *mesh << flush;
232  }
233 
234  // Clean up and exit
235  delete mesh;
236  return 0;
237 }
238 
239 void trans(const Vector &x, Vector &p)
240 {
241  int nnode = (el_type_ == Element::WEDGE)? 3:4;
242 
243  double phi = 2.0 * M_PI * x[2] / nphi_;
244  double theta = theta0_ + phi * ns_ / nnode;
245 
246  double u = (1.5 * (x[0] + x[1]) - 1.0) * r_;
247  double v = sqrt(0.75) * (x[0] - x[1]) * r_;
248 
249  if (el_type_ == Element::WEDGE)
250  {
251  u = (1.5 * (x[0] + x[1]) - 1.0) * r_;
252  v = sqrt(0.75) * (x[0] - x[1]) * r_;
253  }
254  else
255  {
256  u = M_SQRT2 * (x[1] - 0.5) * r_;
257  v = M_SQRT2 * (x[0] - 0.5) * r_;
258  }
259 
260  p[0] = ( R_ + u * cos(theta) + v * sin(theta)) * cos(phi);
261  p[1] = ( R_ + u * cos(theta) + v * sin(theta)) * sin(phi);
262  p[2] = v * cos(theta) - u * sin(theta);
263 }
void AddHex(const int *vi, int attr=1)
Definition: mesh.cpp:1204
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1188
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:696
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:9864
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
int main(int argc, char *argv[])
Definition: ex1.cpp:62
void RemoveInternalBoundaries()
Definition: mesh.cpp:10020
void AddVertex(const double *)
Definition: mesh.cpp:1155
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3924
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
const Element * GetElement(int i) const
Definition: mesh.hpp:775
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void trans(const Vector &x, Vector &p)
Definition: toroid.cpp:239
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:9913
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2242
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)
Definition: optparser.hpp:76
void AddWedge(const int *vi, int attr=1)
Definition: mesh.cpp:1199
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:690
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual int GetNVertices() const =0
Vector data type.
Definition: vector.hpp:48
void pts(int iphi, int t, double x[])
Abstract data type element.
Definition: element.hpp:28
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:779
bool Good() const
Definition: optparser.hpp:125