-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
In progressive mesh generation, change tag::inherent_orientation into tag::orientation, tag::inherent etc. Introduce constructor Mesh ( tag::import, tag::msh, "filename.msh"). Rename export_msh to export_to_file. Switch to gmsh version 4 format. Introduce constructor Mesh ( tag::cube, ... ). Gather source files under src directory.
- Loading branch information
1 parent
2583125
commit b91e706
Showing
33 changed files
with
3,349 additions
and
3,218 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,125 +1,91 @@ | ||
|
||
// example presented in paragraph 2.15 of the manual | ||
// http://manifem.rd.ciencias.ulisboa.pt/manual-manifem.pdf | ||
// a "bumpy" hemisphere with smooth base (submanifold of a manifold, implicit) | ||
// a 3D star-like shape | ||
|
||
#include "maniFEM.h" | ||
#include "math.h" | ||
|
||
using namespace maniFEM; | ||
using namespace std; | ||
|
||
|
||
int main () | ||
|
||
{ Manifold RR3 ( tag::Euclid, tag::of_dim, 3 ); | ||
Function xyz = RR3 .build_coordinate_system ( tag::Lagrange, tag::of_degree, 1 ); | ||
Function x = xyz [0], y = xyz [1], z = xyz [2]; | ||
|
||
Manifold nut = RR3 .implicit ( x*x + y*y + z*z + 1.5*x*y*z == 1. ); | ||
|
||
int n = 10; | ||
|
||
// first, build the base (a closed curve) | ||
|
||
Manifold base = nut .implicit ( x*x + 3.*z == 0. ); | ||
|
||
Cell S ( tag::vertex ); x (S) = 0.; y (S) = -1.; z (S) = 0.; | ||
Cell E ( tag::vertex ); x (E) = 1.; y (E) = 0.; z (E) = 0.; | ||
Cell N ( tag::vertex ); x (N) = 0.; y (N) = 1.; z (N) = 0.; | ||
Cell W ( tag::vertex ); x (W) = -1.; y (W) = 0.; z (W) = 0.; | ||
// no need to project S and N, they are already on 'base' | ||
base .project (E); base .project (W); | ||
Cell mSW ( tag::vertex ); x (mSW) = -1.; y (mSW) = -1.; z (mSW) = 0.; | ||
Cell mSE ( tag::vertex ); x (mSE) = 1.; y (mSE) = -1.; z (mSE) = 0.; | ||
Cell mNE ( tag::vertex ); x (mNE) = 1.; y (mNE) = 1.; z (mNE) = 0.; | ||
Cell mNW ( tag::vertex ); x (mNW) = -1.; y (mNW) = 1.; z (mNW) = 0.; | ||
base .project ( mSW ); base .project ( mSE ); base .project ( mNE ); base .project ( mNW ); | ||
|
||
// now build eight segments, forming the base | ||
Mesh S_mSW ( tag::segment, S .reverse(), mSW, tag::divided_in, n ); | ||
Mesh S_mSE ( tag::segment, S .reverse(), mSE, tag::divided_in, n ); | ||
Mesh E_mSE ( tag::segment, E .reverse(), mSE, tag::divided_in, n ); | ||
Mesh E_mNE ( tag::segment, E .reverse(), mNE, tag::divided_in, n ); | ||
Mesh N_mNE ( tag::segment, N .reverse(), mNE, tag::divided_in, n ); | ||
Mesh N_mNW ( tag::segment, N .reverse(), mNW, tag::divided_in, n ); | ||
Mesh W_mSW ( tag::segment, W .reverse(), mSW, tag::divided_in, n ); | ||
Mesh W_mNW ( tag::segment, W .reverse(), mNW, tag::divided_in, n ); | ||
|
||
// we are done with the bas, now switch back to 'nut' | ||
nut .set_as_working_manifold(); | ||
|
||
// more points : | ||
Cell up ( tag::vertex ); x (up) = 0.; y (up) = 0.; z (up) = 1.; | ||
// no need to project 'up', it is already on 'nut' | ||
Cell mSup ( tag::vertex ); x (mSup) = 0.; y (mSup) = -1.; z (mSup) = 1.; | ||
Cell mEup ( tag::vertex ); x (mEup) = 1.; y (mEup) = 0.; z (mEup) = 1.; | ||
Cell mNup ( tag::vertex ); x (mNup) = 0.; y (mNup) = 1.; z (mNup) = 1.; | ||
Cell mWup ( tag::vertex ); x (mWup) = -1.; y (mWup) = 0.; z (mWup) = 1.; | ||
nut .project ( mSup ); nut .project ( mEup ); nut .project ( mNup ); nut .project ( mWup ); | ||
Cell mSWup ( tag::vertex ); x (mSWup) = -1.; y (mSWup) = -1.; z (mSWup) = 1.; | ||
Cell mSEup ( tag::vertex ); x (mSEup) = 1.; y (mSEup) = -1.; z (mSEup) = 1.; | ||
Cell mNEup ( tag::vertex ); x (mNEup) = 1.; y (mNEup) = 1.; z (mNEup) = 1.; | ||
Cell mNWup ( tag::vertex ); x (mNWup) = -1.; y (mNWup) = 1.; z (mNWup) = 1.; | ||
nut .project ( mSWup ); nut .project ( mSEup ); | ||
nut .project ( mNEup ); nut .project ( mNWup ); | ||
|
||
// more segments : | ||
Mesh S_mSup ( tag::segment, S .reverse(), mSup, tag::divided_in, n ); | ||
Mesh N_mNup ( tag::segment, N .reverse(), mNup, tag::divided_in, n ); | ||
Mesh E_mEup ( tag::segment, E .reverse(), mEup, tag::divided_in, n ); | ||
Mesh W_mWup ( tag::segment, W .reverse(), mWup, tag::divided_in, n ); | ||
Mesh mSW_mSWup ( tag::segment, mSW .reverse(), mSWup, tag::divided_in, n ); | ||
Mesh mSE_mSEup ( tag::segment, mSE .reverse(), mSEup, tag::divided_in, n ); | ||
Mesh mNE_mNEup ( tag::segment, mNE .reverse(), mNEup, tag::divided_in, n ); | ||
Mesh mNW_mNWup ( tag::segment, mNW .reverse(), mNWup, tag::divided_in, n ); | ||
Mesh up_mSup ( tag::segment, up .reverse(), mSup, tag::divided_in, n ); | ||
Mesh up_mEup ( tag::segment, up .reverse(), mEup, tag::divided_in, n ); | ||
Mesh up_mNup ( tag::segment, up .reverse(), mNup, tag::divided_in, n ); | ||
Mesh up_mWup ( tag::segment, up .reverse(), mWup, tag::divided_in, n ); | ||
Mesh mSup_mSEup ( tag::segment, mSup .reverse(), mSEup, tag::divided_in, n ); | ||
Mesh mSup_mSWup ( tag::segment, mSup .reverse(), mSWup, tag::divided_in, n ); | ||
Mesh mEup_mSEup ( tag::segment, mEup .reverse(), mSEup, tag::divided_in, n ); | ||
Mesh mEup_mNEup ( tag::segment, mEup .reverse(), mNEup, tag::divided_in, n ); | ||
Mesh mNup_mNEup ( tag::segment, mNup .reverse(), mNEup, tag::divided_in, n ); | ||
Mesh mNup_mNWup ( tag::segment, mNup .reverse(), mNWup, tag::divided_in, n ); | ||
Mesh mWup_mSWup ( tag::segment, mWup .reverse(), mSWup, tag::divided_in, n ); | ||
Mesh mWup_mNWup ( tag::segment, mWup .reverse(), mNWup, tag::divided_in, n ); | ||
Cell A ( tag::vertex ); x (A) = -1.; y (A) = -1.; z (A) = -1.; | ||
Cell B ( tag::vertex ); x (B) = 1.; y (B) = -1.; z (B) = -1.; | ||
Cell C ( tag::vertex ); x (C) = 1.; y (C) = -1.; z (C) = 1.; | ||
Cell D ( tag::vertex ); x (D) = -1.; y (D) = -1.; z (D) = 1.; | ||
Cell E ( tag::vertex ); x (E) = -1.; y (E) = 1.; z (E) = -1.; | ||
Cell F ( tag::vertex ); x (F) = 1.; y (F) = 1.; z (F) = -1.; | ||
Cell G ( tag::vertex ); x (G) = 1.; y (G) = 1.; z (G) = 1.; | ||
Cell H ( tag::vertex ); x (H) = -1.; y (H) = 1.; z (H) = 1.; | ||
|
||
const double p = 3.3; // recommended values p > 3. | ||
const double q = (1.-p) / 4.; | ||
|
||
// one can define the same curve by different sets of equations | ||
|
||
// for instance, below we use the equation already introduced in paragraph 2.12 | ||
// and add the equation of a plane | ||
|
||
RR3 .implicit ( x*x + q*(y+z)*(y+z) == 2.-p, y == z ); | ||
|
||
// but this is equivalent to (prove it mathematically, check it out in the computer) : | ||
// RR3 .implicit ( x*x + (1.-p)*y*y == 2.-p, y == z ); | ||
// or : RR3 .implicit ( x*x + (1.-p)*z*z == 2.-p, y == z ); | ||
// or : RR3 .implicit ( x*x + y*y - p*z*z == 2.-p, y == z ); | ||
// or : RR3 .implicit ( x*x - p*y*y + z*z == 2.-p, y == z ); | ||
|
||
Mesh AB ( tag::segment, A .reverse(), B, tag::divided_in, 15 ); | ||
Mesh GH ( tag::segment, G .reverse(), H, tag::divided_in, 15 ); | ||
|
||
RR3 .implicit ( x*x + q*(y-z)*(y-z) == 2.-p, y == -z ); | ||
// or : RR3 .implicit ( x*x + (1.-p)*y*y == 2.-p, y == -z ); | ||
// or : RR3 .implicit ( x*x + (1.-p)*z*z == 2.-p, y == -z ); | ||
// or : RR3 .implicit ( x*x + y*y - p*z*z == 2.-p, y == -z ); | ||
// or : RR3 .implicit ( x*x - p*y*y + z*z == 2.-p, y == -z ); | ||
Mesh CD ( tag::segment, C .reverse(), D, tag::divided_in, 15 ); | ||
Mesh EF ( tag::segment, E .reverse(), F, tag::divided_in, 15 ); | ||
|
||
// below is still another set of equations defining the same kind of one-dimensional manifold | ||
|
||
// now the twelve rectangles : | ||
Mesh rect_S_SE ( tag::rectangle, | ||
S_mSE, mSE_mSEup, mSup_mSEup .reverse(), S_mSup .reverse() ); | ||
Mesh rect_S_SW ( tag::rectangle, | ||
S_mSup, mSup_mSWup, mSW_mSWup .reverse(), S_mSW .reverse() ); | ||
Mesh rect_E_SE ( tag::rectangle, | ||
E_mEup, mEup_mSEup, mSE_mSEup .reverse(), E_mSE .reverse() ); | ||
Mesh rect_E_NE ( tag::rectangle, | ||
E_mNE, mNE_mNEup, mEup_mNEup .reverse(), E_mEup .reverse() ); | ||
Mesh rect_N_NE ( tag::rectangle, | ||
N_mNup, mNup_mNEup, mNE_mNEup .reverse(), N_mNE .reverse() ); | ||
Mesh rect_N_NW ( tag::rectangle, | ||
N_mNW, mNW_mNWup, mNup_mNWup .reverse(), N_mNup .reverse() ); | ||
// yes, this sounds a lot like Carry Grant ... | ||
Mesh rect_W_SW ( tag::rectangle, | ||
W_mSW, mSW_mSWup, mWup_mSWup .reverse(), W_mWup .reverse() ); | ||
Mesh rect_W_NW ( tag::rectangle, | ||
W_mWup, mWup_mNWup, mNW_mNWup .reverse(), W_mNW .reverse() ); | ||
Mesh rect_up_SW ( tag::rectangle, | ||
up_mWup, mWup_mSWup, mSup_mSWup .reverse(), up_mSup .reverse() ); | ||
Mesh rect_up_SE ( tag::rectangle, | ||
up_mSup, mSup_mSEup, mEup_mSEup .reverse(), up_mEup .reverse() ); | ||
Mesh rect_up_NE ( tag::rectangle, | ||
up_mEup, mEup_mNEup, mNup_mNEup .reverse(), up_mNup .reverse() ); | ||
Mesh rect_up_NW ( tag::rectangle, | ||
up_mNup, mNup_mNWup, mWup_mNWup .reverse(), up_mWup .reverse() ); | ||
|
||
// and finally join the rectangles : | ||
Mesh bumpy ( tag::join, | ||
{ rect_S_SE, rect_S_SW, rect_E_SE, rect_E_NE, rect_N_NE, rect_N_NW, | ||
rect_W_SW, rect_W_NW, rect_up_SW, rect_up_SE, rect_up_NE, rect_up_NW } ); | ||
RR3 .implicit ( x*x + y*y - p*z*z == 2.-p, - p*x*x + y*y + z*z == 2.-p ); | ||
Mesh AE ( tag::segment, A .reverse(), E, tag::divided_in, 15 ); | ||
Mesh BF ( tag::segment, B .reverse(), F, tag::divided_in, 15 ); | ||
Mesh CG ( tag::segment, C .reverse(), G, tag::divided_in, 15 ); | ||
Mesh DH ( tag::segment, D .reverse(), H, tag::divided_in, 15 ); | ||
|
||
RR3 .implicit ( - p*x*x + y*y + z*z == 2.-p, x*x - p*y*y + z*z == 2.-p ); | ||
Mesh BC ( tag::segment, B .reverse(), C, tag::divided_in, 15 ); | ||
Mesh DA ( tag::segment, D .reverse(), A, tag::divided_in, 15 ); | ||
Mesh FG ( tag::segment, F .reverse(), G, tag::divided_in, 15 ); | ||
Mesh HE ( tag::segment, H .reverse(), E, tag::divided_in, 15 ); | ||
|
||
// for surfaces, we do not have so much freedom : | ||
|
||
bumpy .export_to_file ( tag::msh, "bumpy.msh"); | ||
RR3 .implicit ( x*x + y*y - p*z*z == 2.-p ); | ||
Mesh AEFB ( tag::rectangle, AE, EF, BF .reverse(), AB .reverse() ); | ||
Mesh CGHD ( tag::rectangle, CG, GH, DH .reverse(), CD .reverse() ); | ||
|
||
cout << "produced file bumpy.msh" << endl; | ||
RR3 .implicit ( x*x - p*y*y + z*z == 2.-p ); | ||
Mesh ABCD ( tag::rectangle, AB, BC, CD, DA ); | ||
Mesh EFGH ( tag::rectangle, EF, FG, GH, HE ); | ||
|
||
RR3 .implicit ( - p*x*x + y*y + z*z == 2.-p ); | ||
Mesh DHEA ( tag::rectangle, DH, HE, AE .reverse(), DA .reverse() ); | ||
Mesh BFGC ( tag::rectangle, BF, FG, CG .reverse(), BC .reverse() ); | ||
|
||
// back to the initial 3D Euclidian space : | ||
|
||
RR3 .set_as_working_manifold(); | ||
Mesh cube ( tag::cube, ABCD, EFGH .reverse(), BFGC, DHEA, CGHD, AEFB ); | ||
|
||
cube .export_to_file ( tag::msh, "cube.msh"); | ||
|
||
std::cout << "produced file cube.msh" << std::endl; | ||
|
||
} // end of main | ||
|
||
|
Oops, something went wrong.