MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
klein-bottle.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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// 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
32using namespace std;
33using namespace mfem;
34
35void figure8_trans(const Vector &x, Vector &p);
36void bottle_trans(const Vector &x, Vector &p);
37void bottle2_trans(const Vector &x, Vector &p);
38
39int 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 // The mesh could use quads (default) or triangles
76 // Element::Type el_type = Element::TRIANGLE;
77 Mesh mesh = Mesh::MakeCartesian2D(nx, ny, el_type, 1, 2*M_PI, 2*M_PI);
78
79 mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
80
81 {
82 Array<int> v2v(mesh.GetNV());
83 for (int i = 0; i < v2v.Size(); i++)
84 {
85 v2v[i] = i;
86 }
87 // identify vertices on horizontal lines (without a flip)
88 for (int i = 0; i <= nx; i++)
89 {
90 int v_old = i + ny * (nx + 1);
91 int v_new = i;
92 v2v[v_old] = v_new;
93 }
94 // identify vertices on vertical lines (with a flip)
95 for (int j = 0; j <= ny; j++)
96 {
97 int v_old = nx + j * (nx + 1);
98 int v_new = (ny - j) * (nx + 1);
99 v2v[v_old] = v2v[v_new];
100 }
101 // renumber elements
102 for (int i = 0; i < mesh.GetNE(); i++)
103 {
104 Element *el = mesh.GetElement(i);
105 int *v = el->GetVertices();
106 int nv = el->GetNVertices();
107 for (int j = 0; j < nv; j++)
108 {
109 v[j] = v2v[v[j]];
110 }
111 }
112 // renumber boundary elements
113 for (int i = 0; i < mesh.GetNBE(); i++)
114 {
115 Element *el = mesh.GetBdrElement(i);
116 int *v = el->GetVertices();
117 int nv = el->GetNVertices();
118 for (int j = 0; j < nv; j++)
119 {
120 v[j] = v2v[v[j]];
121 }
122 }
125 }
126
127 switch (trans_type)
128 {
129 case 0: mesh.Transform(figure8_trans); break;
130 case 1: mesh.Transform(bottle_trans); break;
131 case 2: mesh.Transform(bottle2_trans); break;
132 default: mesh.Transform(bottle_trans); break;
133 }
134
135 if (!dg_mesh)
136 {
137 mesh.SetCurvature(order, false, 3, Ordering::byVDIM);
138 }
139
140 GridFunction &nodes = *mesh.GetNodes();
141 for (int i = 0; i < nodes.Size(); i++)
142 {
143 if (std::abs(nodes(i)) < 1e-12)
144 {
145 nodes(i) = 0.0;
146 }
147 }
148
149 ofstream ofs(new_mesh_file);
150 ofs.precision(8);
151 mesh.Print(ofs);
152 ofs.close();
153
154 if (visualization)
155 {
156 char vishost[] = "localhost";
157 int visport = 19916;
158 socketstream sol_sock(vishost, visport);
159 sol_sock.precision(8);
160 sol_sock << "mesh\n" << mesh << flush;
161 }
162
163 return 0;
164}
165
166void figure8_trans(const Vector &x, Vector &p)
167{
168 const real_t r = 2.5;
169 real_t a = r + cos(x(0)/2) * sin(x(1)) - sin(x(0)/2) * sin(2*x(1));
170
171 p.SetSize(3);
172 p(0) = a * cos(x(0));
173 p(1) = a * sin(x(0));
174 p(2) = sin(x(0)/2) * sin(x(1)) + cos(x(0)/2) * sin(2*x(1));
175}
176
177void bottle_trans(const Vector &x, Vector &p)
178{
179 real_t u = x(0);
180 real_t v = x(1) + M_PI_2;
181 real_t a = 6.*cos(u)*(1.+sin(u));
182 real_t b = 16.*sin(u);
183 real_t r = 4.*(1.-cos(u)/2.);
184
185 if (u <= M_PI)
186 {
187 p(0) = a+r*cos(u)*cos(v);
188 p(1) = b+r*sin(u)*cos(v);
189 }
190 else
191 {
192 p(0) = a+r*cos(v+M_PI);
193 p(1) = b;
194 }
195 p(2) = r*sin(v);
196}
197
198void bottle2_trans(const Vector &x, Vector &p)
199{
200 real_t u = x(1)-M_PI_2, v = 2*x(0);
201 const real_t pi = M_PI;
202
203 p(0) = (v<pi ? (2.5-1.5*cos(v))*cos(u) :
204 (v<2*pi ? (2.5-1.5*cos(v))*cos(u) :
205 (v<3*pi ? -2+(2+cos(u))*cos(v) : -2+2*cos(v)-cos(u))));
206 p(1) = (v<pi ? (2.5-1.5*cos(v))*sin(u) :
207 (v<2*pi ? (2.5-1.5*cos(v))*sin(u) :
208 (v<3*pi ? sin(u) : sin(u))));
209 p(2) = (v<pi ? -2.5*sin(v) :
210 (v<2*pi ? 3*v-3*pi :
211 (v<3*pi ? (2+cos(u))*sin(v)+3*pi : -3*v+12*pi)));
212}
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
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
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
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8973
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
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
Definition: mesh.cpp:4243
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 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 ...
Definition: optparser.cpp:151
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
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.
Definition: optparser.hpp:159
Vector data type.
Definition: vector.hpp:80
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
int visport
char vishost[]
int main()
void figure8_trans(const Vector &x, Vector &p)
void bottle_trans(const Vector &x, Vector &p)
void bottle2_trans(const Vector &x, Vector &p)
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
real_t u(const Vector &xvec)
Definition: lor_mms.hpp:22
float real_t
Definition: config.hpp:43
real_t p(const Vector &x, real_t t)
Definition: navier_mms.cpp:53