From fbcd55760c55489124282c3fb0257a087d875444 Mon Sep 17 00:00:00 2001
From: Brian Helenbrook <bhelenbr@clarkson.edu>
Date: Tue, 18 Nov 2008 15:05:18 +0000
Subject: [PATCH] Added uniform refinement routines
---
Common/CommandLine.cpp | 5 +
Fltk/Callbacks.cpp | 8 ++
Fltk/Callbacks.h | 1 +
Fltk/GUI.cpp | 5 +-
Mesh/HighOrder.h | 2 +
Mesh/Refine.cpp | 245 +++++++++++++++++++++++++++++++++++++++++
contrib/NR/Makefile | 2 +-
7 files changed, 265 insertions(+), 3 deletions(-)
create mode 100644 Mesh/Refine.cpp
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index a24d686c2f..fd89bc4fd8 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -42,6 +42,7 @@ void Print_Usage(const char *name)
Msg::Direct(" -tol float Set geometrical tolerance");
Msg::Direct("Mesh options:");
Msg::Direct(" -1, -2, -3 Perform 1D, 2D or 3D mesh generation, then exit");
+ Msg::Direct(" -refine Perform uniform mesh refinement, then exit");
Msg::Direct(" -part Partition after batch mesh generation");
Msg::Direct(" -saveall Save all elements (discard physical group definitions)");
Msg::Direct(" -o file Specify mesh output file name");
@@ -202,6 +203,10 @@ void Get_Options(int argc, char *argv[])
else
Msg::Fatal("Missing string");
CTX.batch = 5;
+ }
+ else if(!strcmp(argv[i] + 1, "refine")) {
+ CTX.batch = 6;
+ i++;
}
else if(!strcmp(argv[i] + 1, "part")) {
CTX.batchAfterMesh = 1;
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 40b7c46eea..d9046f8409 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -3826,6 +3826,14 @@ void mesh_optimize_cb(CALLBACK_ARGS)
Msg::StatusBar(2, false, " ");
}
+void mesh_refine_cb(CALLBACK_ARGS)
+{
+ Refine(GModel::current(), CTX.mesh.second_order_linear);
+ CTX.mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
+ Draw();
+ Msg::StatusBar(2, false, " ");
+}
+
void mesh_optimize_netgen_cb(CALLBACK_ARGS)
{
if(CTX.threads_lock) {
diff --git a/Fltk/Callbacks.h b/Fltk/Callbacks.h
index d9a199a05d..f1b9ffad22 100644
--- a/Fltk/Callbacks.h
+++ b/Fltk/Callbacks.h
@@ -232,6 +232,7 @@ void mesh_remesh_cb(CALLBACK_ARGS);
void mesh_update_edges_cb(CALLBACK_ARGS);
void mesh_parameterize_cb(CALLBACK_ARGS);
void mesh_degree_cb(CALLBACK_ARGS);
+void mesh_refine_cb(CALLBACK_ARGS);
void mesh_optimize_cb(CALLBACK_ARGS);
void mesh_optimize_netgen_cb(CALLBACK_ARGS);
void mesh_partition_cb(CALLBACK_ARGS);
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 4c9fb61cb1..78d7e2c264 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -304,6 +304,7 @@ Context_Item menu_mesh[] = {
{"3D", (Fl_Callback *)mesh_3d_cb} ,
{"First order", (Fl_Callback *)mesh_degree_cb, (void*)1 } ,
{"Second order", (Fl_Callback *)mesh_degree_cb, (void*)2 } ,
+ {"Refine", (Fl_Callback *)mesh_refine_cb} ,
{"Optimize", (Fl_Callback *)mesh_optimize_cb} ,
#if defined(HAVE_NETGEN)
{"Optimize (Netgen)", (Fl_Callback *)mesh_optimize_netgen_cb} ,
@@ -4819,8 +4820,8 @@ void GUI::create_about_window()
sprintf(buffer, "@c@.Build options: %s", lines[i].c_str());
else
sprintf(buffer, "@c@.%s", lines[i].c_str());
- o->add(buffer);
- }
+ o->add(buffer);
+ }
sprintf(buffer, "@c@.Packaged by: %s", Get_GmshPackager());
o->add(buffer);
o->add(" ");
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index 98d2d62653..6e1d1a678d 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -13,4 +13,6 @@ void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false);
void checkHighOrderTriangles(GModel *m);
bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m);
+void Refine(GModel *m, bool linear=true);
+
#endif
diff --git a/Mesh/Refine.cpp b/Mesh/Refine.cpp
new file mode 100644
index 0000000000..ccefea7346
--- /dev/null
+++ b/Mesh/Refine.cpp
@@ -0,0 +1,245 @@
+/*
+ * refine.cpp
+ * GMSH
+ *
+ * Created by Brian Helenbrook on 11/13/08.
+ * Copyright 2008 Clarkson University. All rights reserved.
+ *
+ */
+#include "HighOrder.h"
+#include "meshGFaceOptimize.h"
+#include "MElement.h"
+#include "GmshMessage.h"
+#include "OS.h"
+#include "Numeric.h"
+#include "Context.h"
+#include "GmshMatrix.h"
+#include "FunctionSpace.h"
+
+void Subdivide(GEdge *ge)
+{
+ std::vector<MLine *> lines2;
+ for(unsigned int i = 0; i < ge->lines.size(); i++){
+ MLine3 *l = dynamic_cast<MLine3 *>(ge->lines[i]);
+ lines2.push_back(new MLine(l->getVertex(0), l->getVertex(2)));
+ lines2.push_back(new MLine(l->getVertex(2), l->getVertex(1)));
+ delete l;
+ }
+ ge->lines = lines2;
+ ge->deleteVertexArrays();
+}
+
+void Subdivide(GFace *gf)
+{
+ std::vector<MTriangle*> triangles2;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle6 *t = dynamic_cast<MTriangle6 *>(gf->triangles[i]);
+ triangles2.push_back
+ (new MTriangle(t->getVertex(0), t->getVertex(3), t->getVertex(5)));
+ triangles2.push_back
+ (new MTriangle(t->getVertex(3), t->getVertex(4), t->getVertex(5)));
+ triangles2.push_back
+ (new MTriangle(t->getVertex(3), t->getVertex(1), t->getVertex(4)));
+ triangles2.push_back
+ (new MTriangle(t->getVertex(5), t->getVertex(4), t->getVertex(2)));
+ delete t;
+ }
+ gf->triangles = triangles2;
+
+ std::vector<MQuadrangle*> quadrangles2;
+ for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
+ MQuadrangle9 *q = dynamic_cast<MQuadrangle9 *>(gf->quadrangles[i]);
+ quadrangles2.push_back
+ (new MQuadrangle(q->getVertex(0), q->getVertex(4), q->getVertex(8), q->getVertex(7)));
+ quadrangles2.push_back
+ (new MQuadrangle(q->getVertex(4), q->getVertex(1), q->getVertex(5), q->getVertex(8)));
+ quadrangles2.push_back
+ (new MQuadrangle(q->getVertex(8), q->getVertex(5), q->getVertex(2), q->getVertex(6)));
+ quadrangles2.push_back
+ (new MQuadrangle(q->getVertex(7), q->getVertex(8), q->getVertex(6), q->getVertex(3)));
+
+ delete q;
+ }
+ gf->quadrangles = quadrangles2;
+ gf->deleteVertexArrays();
+}
+
+void Subdivide(GRegion *gr)
+{
+
+ std::vector<MTetrahedron*> tetrahedra2;
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ MTetrahedron10 *t = dynamic_cast<MTetrahedron10 *>(gr->tetrahedra[i]);
+
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(0), t->getVertex(4), t->getVertex(7), t->getVertex(6)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(1), t->getVertex(4), t->getVertex(5), t->getVertex(9)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(2), t->getVertex(5), t->getVertex(6), t->getVertex(8)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(3), t->getVertex(7), t->getVertex(9), t->getVertex(8)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(5), t->getVertex(8), t->getVertex(7), t->getVertex(9)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(5), t->getVertex(7), t->getVertex(4), t->getVertex(9)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(7), t->getVertex(8), t->getVertex(5), t->getVertex(6)));
+ tetrahedra2.push_back
+ (new MTetrahedron(t->getVertex(4), t->getVertex(7), t->getVertex(5), t->getVertex(6)));
+ delete t;
+ }
+ gr->tetrahedra = tetrahedra2;
+
+ std::vector<MHexahedron*> hexahedra2;
+ for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
+ MHexahedron27 *h = dynamic_cast<MHexahedron27 *>(gr->hexahedra[i]);
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(0), h->getVertex(8), h->getVertex(20), h->getVertex(9),
+ h->getVertex(10), h->getVertex(21), h->getVertex(26), h->getVertex(22)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(10), h->getVertex(21), h->getVertex(26), h->getVertex(22),
+ h->getVertex(4), h->getVertex(16), h->getVertex(25), h->getVertex(17)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(8), h->getVertex(1), h->getVertex(11), h->getVertex(20),
+ h->getVertex(21), h->getVertex(12), h->getVertex(23), h->getVertex(26)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(21), h->getVertex(12), h->getVertex(23), h->getVertex(26),
+ h->getVertex(16), h->getVertex(5), h->getVertex(18), h->getVertex(25)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(9), h->getVertex(20), h->getVertex(13), h->getVertex(3),
+ h->getVertex(22), h->getVertex(26), h->getVertex(24), h->getVertex(15)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(22), h->getVertex(26), h->getVertex(24), h->getVertex(15),
+ h->getVertex(17), h->getVertex(25), h->getVertex(19), h->getVertex(7)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(20), h->getVertex(11), h->getVertex(2), h->getVertex(13),
+ h->getVertex(26), h->getVertex(23), h->getVertex(14), h->getVertex(24)));
+ hexahedra2.push_back
+ (new MHexahedron(h->getVertex(26), h->getVertex(23), h->getVertex(14), h->getVertex(24),
+ h->getVertex(25), h->getVertex(18), h->getVertex(6), h->getVertex(19)));
+ delete h;
+ }
+ gr->hexahedra = hexahedra2;
+
+ std::vector<MPrism*> prisms2;
+ for(unsigned int i = 0; i < gr->prisms.size(); i++){
+ MPrism18 *p = dynamic_cast<MPrism18 *>(gr->prisms[i]);
+
+ prisms2.push_back
+ (new MPrism(p->getVertex(0), p->getVertex(6), p->getVertex(7),
+ p->getVertex(8), p->getVertex(15), p->getVertex(16)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(8), p->getVertex(15), p->getVertex(16),
+ p->getVertex(3), p->getVertex(12), p->getVertex(13)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(6), p->getVertex(1), p->getVertex(9),
+ p->getVertex(15), p->getVertex(10), p->getVertex(17)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(15), p->getVertex(10), p->getVertex(17),
+ p->getVertex(12), p->getVertex(4), p->getVertex(14)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(7), p->getVertex(9), p->getVertex(2),
+ p->getVertex(16), p->getVertex(17), p->getVertex(11)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(16), p->getVertex(17), p->getVertex(11),
+ p->getVertex(13), p->getVertex(14), p->getVertex(5)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(9), p->getVertex(7), p->getVertex(6),
+ p->getVertex(17), p->getVertex(16), p->getVertex(15)));
+ prisms2.push_back
+ (new MPrism(p->getVertex(17), p->getVertex(16), p->getVertex(15),
+ p->getVertex(14), p->getVertex(13), p->getVertex(12)));
+
+ delete p;
+ }
+ gr->prisms = prisms2;
+
+
+ std::vector<MPyramid*> pyramids2;
+ for(unsigned int i = 0; i < gr->pyramids.size(); i++){
+ MPyramid14 *p = dynamic_cast<MPyramid14 *>(gr->pyramids[i]);
+
+ /* BASE */
+ pyramids2.push_back
+ (new MPyramid(p->getVertex(0), p->getVertex(5), p->getVertex(13), p->getVertex(6), p->getVertex(7)));
+ pyramids2.push_back
+ (new MPyramid(p->getVertex(5), p->getVertex(1), p->getVertex(8), p->getVertex(13), p->getVertex(9)));
+ pyramids2.push_back
+ (new MPyramid(p->getVertex(13), p->getVertex(8), p->getVertex(2), p->getVertex(10), p->getVertex(11)));
+ pyramids2.push_back
+ (new MPyramid(p->getVertex(6), p->getVertex(13), p->getVertex(10), p->getVertex(3), p->getVertex(12)));
+
+ /* Split remaining into tets */
+ /* Top */
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(7), p->getVertex(9), p->getVertex(12), p->getVertex(4))));
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(9), p->getVertex(11), p->getVertex(12), p->getVertex(4))));
+
+ /* Upside down one */
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(9), p->getVertex(12), p->getVertex(11), p->getVertex(13))));
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(7), p->getVertex(12), p->getVertex(9), p->getVertex(13))));
+
+
+ /* Four tets around bottom perimeter */
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(7), p->getVertex(9), p->getVertex(5), p->getVertex(13))));
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(9), p->getVertex(11), p->getVertex(8), p->getVertex(13))));
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(12), p->getVertex(10), p->getVertex(11), p->getVertex(13))));
+ gr->tetrahedra.push_back(
+ (new MTetrahedron(p->getVertex(7), p->getVertex(6), p->getVertex(12), p->getVertex(13))));
+ delete p;
+ }
+ gr->pyramids = pyramids2;
+ gr->deleteVertexArrays();
+}
+
+
+
+
+
+
+void Refine(GModel *m, bool linear)
+{
+ // refine all edges in a mesh by inserting a point on the midpoint
+ // then recreate edge, face, and region definitions
+ //
+ // - if linear is set to true, new vertices are created by linear
+ // interpolation between existing ones. If not, new vertices are
+ // created on the exact geometry, provided that the geometrical
+ // edges/faces are discretized with 1D/2D elements. (I.e., if
+ // there are only 3D elements in the mesh then any new vertices on
+ // the boundary will always be created by linear interpolation,
+ // whether linear is set or not.)
+
+
+#if !defined(HAVE_GSL)
+ Msg::Error("Mesh refinement requires the GSL");
+ return;
+#endif
+
+ Msg::StatusBar(1, true, "Generating refined mesh ...");
+ double t1 = Cpu();
+
+ // first, set order to 2 to generate vertex positions
+ SetOrderN(m, 2, linear, false);
+
+ // Now subdivide entities to create linear mesh
+ for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
+ Subdivide(*it);
+
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ Subdivide(*it);
+
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
+ Subdivide(*it);
+
+ double t2 = Cpu();
+ Msg::Info("Mesh refinement complete (%g s)", t2 - t1);
+ Msg::StatusBar(1, true, "Mesh");
+}
diff --git a/contrib/NR/Makefile b/contrib/NR/Makefile
index a3e8bcad58..5438a4f1fd 100644
--- a/contrib/NR/Makefile
+++ b/contrib/NR/Makefile
@@ -58,4 +58,4 @@ lnsrch.o: lnsrch.cpp nrutil.h
lubksb.o: lubksb.cpp
ludcmp.o: ludcmp.cpp nrutil.h
newt.o: newt.cpp nrutil.h
-nrutil.o: nrutil.cpp ../../Common/Message.h
+nrutil.o: nrutil.cpp ../../Common/GmshMessage.h
--
GitLab