MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mobius-strip.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// Mobius Strip Miniapp: Generate Mobius strip meshes
14// ---------------------------------------------------
15//
16// This miniapp generates various Mobius strip-like surface meshes. It is a good
17// way to generate complex surface meshes. Manipulating the mesh topology and
18// performing mesh transformation are demonstrated. The mobius-strip mesh in the
19// data/ directory was generated with this miniapp.
20//
21// Compile with: make mobius-strip
22//
23// Sample runs: mobius-strip
24// mobius-strip -t 4.5 -nx 16
25// mobius-strip -c 1 -t 1
26// mobius-strip -c 1 -t 4 -nx 16
27// mobius-strip -c 0 -t 0.75
28
29#include "mfem.hpp"
30#include <fstream>
31#include <iostream>
32
33using namespace std;
34using namespace mfem;
35
37void mobius_trans(const Vector &x, Vector &p);
38
39int main(int argc, char *argv[])
40{
41 const char *new_mesh_file = "mobius-strip.mesh";
42 int nx = 8;
43 int ny = 2;
44 int order = 3;
45 int close_strip = 2;
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(&close_strip, "-c", "--close-strip",
59 "How to close the strip: 0 - open, 1 - closed, 2 - twisted.");
60 args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
61 "Use discontinuous or continuous space for the mesh nodes.");
62 args.AddOption(&num_twists, "-t", "--num-twists",
63 "Number of twists of the strip.");
64 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
65 "--no-visualization",
66 "Enable or disable GLVis visualization.");
67 args.Parse();
68 if (!args.Good())
69 {
70 args.PrintUsage(cout);
71 return 1;
72 }
73 args.PrintOptions(cout);
74
75 // The mesh could use quads (default) or triangles
77 // Element::Type el_type = Element::TRIANGLE;
78 Mesh mesh = Mesh::MakeCartesian2D(nx, ny, el_type, 1, 2*M_PI, 2.0);
79
80 mesh.SetCurvature(order, true, 3, Ordering::byVDIM);
81
82 if (close_strip)
83 {
84 Array<int> v2v(mesh.GetNV());
85 for (int i = 0; i < v2v.Size(); i++)
86 {
87 v2v[i] = i;
88 }
89 // identify vertices on vertical lines (with a flip)
90 for (int j = 0; j <= ny; j++)
91 {
92 int v_old = nx + j * (nx + 1);
93 int v_new = ((close_strip == 1) ? j : (ny - j)) * (nx + 1);
94 v2v[v_old] = v_new;
95 }
96 // renumber elements
97 for (int i = 0; i < mesh.GetNE(); i++)
98 {
99 Element *el = mesh.GetElement(i);
100 int *v = el->GetVertices();
101 int nv = el->GetNVertices();
102 for (int j = 0; j < nv; j++)
103 {
104 v[j] = v2v[v[j]];
105 }
106 }
107 // renumber boundary elements
108 for (int i = 0; i < mesh.GetNBE(); i++)
109 {
110 Element *el = mesh.GetBdrElement(i);
111 int *v = el->GetVertices();
112 int nv = el->GetNVertices();
113 for (int j = 0; j < nv; j++)
114 {
115 v[j] = v2v[v[j]];
116 }
117 }
120 }
121
123
124 if (!dg_mesh)
125 {
126 mesh.SetCurvature(order, false, 3, Ordering::byVDIM);
127 }
128
129 GridFunction &nodes = *mesh.GetNodes();
130 for (int i = 0; i < nodes.Size(); i++)
131 {
132 if (std::abs(nodes(i)) < 1e-12)
133 {
134 nodes(i) = 0.0;
135 }
136 }
137
138 ofstream ofs(new_mesh_file);
139 ofs.precision(8);
140 mesh.Print(ofs);
141 ofs.close();
142
143 if (visualization)
144 {
145 char vishost[] = "localhost";
146 int visport = 19916;
147 socketstream sol_sock(vishost, visport);
148 sol_sock.precision(8);
149 sol_sock << "mesh\n" << mesh << flush;
150 }
151
152 return 0;
153}
154
155void mobius_trans(const Vector &x, Vector &p)
156{
157 real_t a = 1.0 + 0.5 * (x[1] - 1.0) * cos( num_twists * x[0] );
158
159 p.SetSize(3);
160 p[0] = a * cos( x[0] );
161 p[1] = a * sin( x[0] );
162 p[2] = 0.5 * (x[1] - 1.0) * sin( num_twists * x[0] );
163}
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()
real_t a
Definition: lissajous.cpp:41
void mobius_trans(const Vector &x, Vector &p)
real_t num_twists
float real_t
Definition: config.hpp:43
real_t p(const Vector &x, real_t t)
Definition: navier_mms.cpp:53