MFEM  v3.3
Finite element discretization library
klein-bottle.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 //
12 // ---------------------------------------------------
13 // Klein Bottle Miniapp: Generate Klein bottle meshes
14 // ---------------------------------------------------
15 //
16 // This miniapp generates three types of Klein bottle surfaces. It is similar to
17 // the mobius-strip miniapp. The klein-bottle and klein-donut meshes in the
18 // data/ directory were generated with this miniapp.
19 //
20 // Compile with: make klein-bottle
21 //
22 // Sample runs: klein-bottle
23 // klein-bottle -o 6 -nx 8 -ny 4
24 // klein-bottle -t 0
25 // klein-bottle -t 0 -o 6 -nx 6 -ny 4
26 // klein-bottle -t 2
27 
28 #include "mfem.hpp"
29 #include <fstream>
30 #include <iostream>
31 
32 using namespace std;
33 using namespace mfem;
34 
35 void figure8_trans(const Vector &x, Vector &p);
36 void bottle_trans(const Vector &x, Vector &p);
37 void bottle2_trans(const Vector &x, Vector &p);
38 
39 int main(int argc, char *argv[])
40 {
41  const char *new_mesh_file = "klein-bottle.mesh";
42  int nx = 16;
43  int ny = 8;
44  int order = 3;
45  int trans_type = 1;
46  bool dg_mesh = false;
47  bool visualization = true;
48 
49  OptionsParser args(argc, argv);
50  args.AddOption(&new_mesh_file, "-m", "--mesh-out-file",
51  "Output Mesh file to write.");
52  args.AddOption(&nx, "-nx", "--num-elements-x",
53  "Number of elements in x-direction.");
54  args.AddOption(&ny, "-ny", "--num-elements-y",
55  "Number of elements in y-direction.");
56  args.AddOption(&order, "-o", "--mesh-order",
57  "Order (polynomial degree) of the mesh elements.");
58  args.AddOption(&trans_type, "-t", "--transformation-type",
59  "Set the transformation type: 0 - \"figure-8\","
60  " 1 - \"bottle\", 2 - \"bottle2\".");
61  args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
62  "Use discontinuous or continuous space for the mesh nodes.");
63  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
64  "--no-visualization",
65  "Enable or disable GLVis visualization.");
66  args.Parse();
67  if (!args.Good())
68  {
69  args.PrintUsage(cout);
70  return 1;
71  }
72  args.PrintOptions(cout);
73 
74  Mesh *mesh;
75  // The mesh could use quads (default) or triangles
76  Element::Type el_type = Element::QUADRILATERAL;
77  // Element::Type el_type = Element::TRIANGLE;
78  mesh = new Mesh(nx, ny, el_type, 1, 2*M_PI, 2*M_PI);
79 
80  mesh->SetCurvature(order, true, 3, Ordering::byVDIM);
81 
82  {
83  Array<int> v2v(mesh->GetNV());
84  for (int i = 0; i < v2v.Size(); i++)
85  {
86  v2v[i] = i;
87  }
88  // identify vertices on horizontal lines (without a flip)
89  for (int i = 0; i <= nx; i++)
90  {
91  int v_old = i + ny * (nx + 1);
92  int v_new = i;
93  v2v[v_old] = v_new;
94  }
95  // identify vertices on vertical lines (with a flip)
96  for (int j = 0; j <= ny; j++)
97  {
98  int v_old = nx + j * (nx + 1);
99  int v_new = (ny - j) * (nx + 1);
100  v2v[v_old] = v2v[v_new];
101  }
102  // renumber elements
103  for (int i = 0; i < mesh->GetNE(); i++)
104  {
105  Element *el = mesh->GetElement(i);
106  int *v = el->GetVertices();
107  int nv = el->GetNVertices();
108  for (int j = 0; j < nv; j++)
109  {
110  v[j] = v2v[v[j]];
111  }
112  }
113  // renumber boundary elements
114  for (int i = 0; i < mesh->GetNBE(); i++)
115  {
116  Element *el = mesh->GetBdrElement(i);
117  int *v = el->GetVertices();
118  int nv = el->GetNVertices();
119  for (int j = 0; j < nv; j++)
120  {
121  v[j] = v2v[v[j]];
122  }
123  }
124  mesh->RemoveUnusedVertices();
125  mesh->RemoveInternalBoundaries();
126  }
127 
128  switch (trans_type)
129  {
130  case 0: mesh->Transform(figure8_trans); break;
131  case 1: mesh->Transform(bottle_trans); break;
132  case 2: mesh->Transform(bottle2_trans); break;
133  default: mesh->Transform(bottle_trans); break;
134  }
135 
136  if (!dg_mesh)
137  {
138  mesh->SetCurvature(order, false, 3, Ordering::byVDIM);
139  }
140 
141  GridFunction &nodes = *mesh->GetNodes();
142  for (int i = 0; i < nodes.Size(); i++)
143  {
144  if (std::abs(nodes(i)) < 1e-12)
145  {
146  nodes(i) = 0.0;
147  }
148  }
149 
150  ofstream ofs(new_mesh_file);
151  ofs.precision(8);
152  mesh->Print(ofs);
153  ofs.close();
154 
155  if (visualization)
156  {
157  char vishost[] = "localhost";
158  int visport = 19916;
159  socketstream sol_sock(vishost, visport);
160  sol_sock.precision(8);
161  sol_sock << "mesh\n" << *mesh << flush;
162  }
163 
164  delete mesh;
165  return 0;
166 }
167 
168 void figure8_trans(const Vector &x, Vector &p)
169 {
170  const double r = 2.5;
171  double a = r + cos(x(0)/2) * sin(x(1)) - sin(x(0)/2) * sin(2*x(1));
172 
173  p.SetSize(3);
174  p(0) = a * cos(x(0));
175  p(1) = a * sin(x(0));
176  p(2) = sin(x(0)/2) * sin(x(1)) + cos(x(0)/2) * sin(2*x(1));
177 }
178 
179 void bottle_trans(const Vector &x, Vector &p)
180 {
181  double u = x(0);
182  double v = x(1) + M_PI_2;
183  double a = 6.*cos(u)*(1.+sin(u));
184  double b = 16.*sin(u);
185  double r = 4.*(1.-cos(u)/2.);
186 
187  if (u <= M_PI)
188  {
189  p(0) = a+r*cos(u)*cos(v);
190  p(1) = b+r*sin(u)*cos(v);
191  }
192  else
193  {
194  p(0) = a+r*cos(v+M_PI);
195  p(1) = b;
196  }
197  p(2) = r*sin(v);
198 }
199 
200 void bottle2_trans(const Vector &x, Vector &p)
201 {
202  double u = x(1)-M_PI_2, v = 2*x(0);
203  const double pi = M_PI;
204 
205  p(0) = (v<pi ? (2.5-1.5*cos(v))*cos(u) :
206  (v<2*pi ? (2.5-1.5*cos(v))*cos(u) :
207  (v<3*pi ? -2+(2+cos(u))*cos(v) : -2+2*cos(v)-cos(u))));
208  p(1) = (v<pi ? (2.5-1.5*cos(v))*sin(u) :
209  (v<2*pi ? (2.5-1.5*cos(v))*sin(u) :
210  (v<3*pi ? sin(u) : sin(u))));
211  p(2) = (v<pi ? -2.5*sin(v) :
212  (v<2*pi ? 3*v-3*pi :
213  (v<3*pi ? (2+cos(u))*sin(v)+3*pi : -3*v+12*pi)));
214 }
void bottle2_trans(const Vector &x, Vector &p)
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:310
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:587
int Size() const
Returns the size of the vector.
Definition: vector.hpp:106
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:8228
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:584
STL namespace.
void RemoveInternalBoundaries()
Definition: mesh.cpp:8384
void figure8_trans(const Vector &x, Vector &p)
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3278
Type
Constants for the classes derived from Element.
Definition: element.hpp:37
const Element * GetElement(int i) const
Definition: mesh.hpp:642
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:8277
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:74
int main(int argc, char *argv[])
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:581
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
virtual int GetNVertices() const =0
Vector data type.
Definition: vector.hpp:36
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5368
void bottle_trans(const Vector &x, Vector &p)
Abstract data type element.
Definition: element.hpp:27
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:646
virtual void Print(std::ostream &out=std::cout) const
Definition: mesh.hpp:988
bool Good() const
Definition: optparser.hpp:120