[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jun 19 21:47:18 CDT 2013
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/a9dbda44c3d3/
Changeset: a9dbda44c3d3
Branch: None
User: iulian07
Date: 2013-06-19 23:00:49
Summary: add MPAS test
it starts with an MPAS initial file, it manufactures a departure MPAS, and computes
intersection between these 2 meshes (that are on the same sphere of radius 1)
it is almost the same as case1_test, except it is an MPAS mesh
Affected #: 4 files
diff --git a/tools/mbcslam/CslamUtils.cpp b/tools/mbcslam/CslamUtils.cpp
index 06c3e18..a819825 100644
--- a/tools/mbcslam/CslamUtils.cpp
+++ b/tools/mbcslam/CslamUtils.cpp
@@ -880,7 +880,7 @@ ErrorCode enforce_convexity(Interface * mb, EntityHandle lset)
// we should create a queue with new polygons that need processing for reflex angles
// (obtuse)
std::queue<EntityHandle> newPolys;
- int brokenQuads=0;
+ int brokenPolys=0;
for (Range::iterator eit = inputRange.begin(); eit != inputRange.end()
|| !newPolys.empty() ; eit++)
{
@@ -974,7 +974,7 @@ ErrorCode enforce_convexity(Interface * mb, EntityHandle lset)
if (MB_SUCCESS != rval)
return rval;
mb->remove_entities(lset, &eh, 1);
- brokenQuads++;
+ brokenPolys++;
/*std::cout<<"remove: " ;
mb->list_entities(&eh, 1);
@@ -985,7 +985,7 @@ ErrorCode enforce_convexity(Interface * mb, EntityHandle lset)
}
}
}
- std::cout << brokenQuads << " quads were decomposed in triangles \n";
+ std::cout << brokenPolys << " concave polygons were decomposed in convex ones \n";
return MB_SUCCESS;
}
diff --git a/tools/mbcslam/Makefile.am b/tools/mbcslam/Makefile.am
index f7785e7..48c122e 100644
--- a/tools/mbcslam/Makefile.am
+++ b/tools/mbcslam/Makefile.am
@@ -39,7 +39,7 @@ libmbcslam_la_includedir = $(includedir)
cfgdir = $(libdir)
TESTS = intx_on_sphere_test intx_in_plane_test spec_visu_test spherical_area_test \
- case1_test
+ case1_test intx_mpas
noinst_PROGRAMS = cslam_par_test
if NETCDF_FILE
noinst_PROGRAMS += process_arm
@@ -51,6 +51,7 @@ intx_in_plane_test_SOURCES = intx_in_plane_test.cpp
spec_visu_test_SOURCES = spec_visu_test.cpp
spherical_area_test_SOURCES = spherical_area_test.cpp
case1_test_SOURCES = case1_test.cpp
+intx_mpas_SOURCES = intx_mpas.cpp
process_arm_SOURCES = process_arm.cpp
cslam_par_test_SOURCES = cslam_par_test.cpp
diff --git a/tools/mbcslam/intx_mpas.cpp b/tools/mbcslam/intx_mpas.cpp
new file mode 100644
index 0000000..f29ff34
--- /dev/null
+++ b/tools/mbcslam/intx_mpas.cpp
@@ -0,0 +1,210 @@
+/*
+ * mpas file test
+ *
+ * Created on: Feb 12, 2013
+ */
+
+// copy from case1 test
+
+#include <iostream>
+#include <sstream>
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include "moab/Core.hpp"
+#include "moab/Interface.hpp"
+#include "Intx2MeshOnSphere.hpp"
+#include <math.h>
+#include "moab/ProgOptions.hpp"
+#include "MBTagConventions.hpp"
+
+#include "CslamUtils.hpp"
+
+// for M_PI
+#include <math.h>
+
+#define STRINGIFY_(X) #X
+#define STRINGIFY(X) STRINGIFY_(X)
+
+using namespace moab;
+// some input data
+double gtol = 1.e-9; // this is for geometry tolerance
+std::string input_mesh_file("mpas.vtk"); // input file, plain vtk, mpas on sphere
+// radius is always 1?
+//double CubeSide = 6.; // the above file starts with cube side 6; radius depends on cube side
+double t = 0.1, delta_t = 0.1; // check the script
+
+ErrorCode manufacture_lagrange_mesh_on_sphere(Interface * mb,
+ EntityHandle euler_set, EntityHandle & lagr_set)
+{
+ ErrorCode rval = MB_SUCCESS;
+
+ /*
+ * get all plys first, then vertices, then move them on the surface of the sphere
+ * radius is 1., always
+ *
+ */
+ double radius = 1.;
+ Range polygons;
+ rval = mb->get_entities_by_dimension(euler_set, 2, polygons);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ Range connecVerts;
+ rval = mb->get_connectivity(polygons, connecVerts);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // create new set
+ rval = mb->create_meshset(MESHSET_SET, lagr_set);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // get the coordinates of the old mesh, and move it around the sphere according to case 1
+ // now put the vertices in the right place....
+ //int vix=0; // vertex index in new array
+
+ // first create departure points (vertices in the lagrange mesh)
+ // then connect them in quads
+ std::map<EntityHandle, EntityHandle> newNodes;
+ for (Range::iterator vit = connecVerts.begin(); vit != connecVerts.end();
+ vit++)
+ {
+ EntityHandle oldV = *vit;
+ CartVect posi;
+ rval = mb->get_coords(&oldV, 1, &(posi[0]));
+ if (MB_SUCCESS != rval)
+ return rval;
+ // cslam utils, case 1
+ CartVect newPos;
+ departure_point_case1(posi, t, delta_t, newPos);
+ newPos = radius * newPos;
+ EntityHandle new_vert;
+ rval = mb->create_vertex(&(newPos[0]), new_vert);
+ if (MB_SUCCESS != rval)
+ return rval;
+ newNodes[oldV] = new_vert;
+ }
+ EntityHandle new_conn[MAXEDGES]; // up to 10, for the time being
+ for (Range::iterator it = polygons.begin(); it != polygons.end(); it++)
+ {
+ EntityHandle q = *it;
+ int nnodes;
+ const EntityHandle * conn4;
+ rval = mb->get_connectivity(q, conn4, nnodes);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ for (int i = 0; i < nnodes; i++)
+ {
+ EntityHandle v1 = conn4[i];
+ new_conn[i] = newNodes[v1];
+ }
+ EntityHandle new_poly;
+ rval = mb->create_element(MBQUAD, new_conn, nnodes, new_poly);
+ if (MB_SUCCESS != rval)
+ return rval;
+ rval = mb->add_entities(lagr_set, &new_poly, 1);
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+
+ return rval;
+}
+int main(int argc, char **argv)
+{
+
+
+ const char *filename_mesh1 = STRINGIFY(SRCDIR) "/mpas.vtk";
+ if (argc > 1)
+ {
+ int index = 1;
+ while (index < argc)
+ {
+ if (!strcmp(argv[index], "-gtol")) // this is for geometry tolerance
+ {
+ gtol = atof(argv[++index]);
+ }
+ if (!strcmp(argv[index], "-dt"))
+ {
+ delta_t = atof(argv[++index]);
+ }
+ if (!strcmp(argv[index], "-input"))
+ {
+ filename_mesh1 = argv[++index];
+ }
+ index++;
+ }
+ }
+ std::cout << " case 1: use -gtol " << gtol << " -dt " << delta_t <<
+ " -input " << filename_mesh1 << "\n";
+
+
+
+ Core moab;
+ Interface & mb = moab;
+ EntityHandle euler_set;
+ ErrorCode rval;
+ rval = mb.create_meshset(MESHSET_SET, euler_set);
+ if (MB_SUCCESS != rval)
+ return 1;
+
+ rval = mb.load_file(filename_mesh1, &euler_set);
+
+ if (MB_SUCCESS != rval)
+ return 1;
+
+ // everybody will get a DP tag, including the non owned entities; so exchange tags is not required for LOC (here)
+ EntityHandle lagrange_set;
+ rval = manufacture_lagrange_mesh_on_sphere(&mb, euler_set, lagrange_set);
+ if (MB_SUCCESS != rval)
+ return 1;
+ rval = mb.write_file("lagrIni.h5m", 0, 0, &lagrange_set, 1);
+ if (MB_SUCCESS != rval)
+ std::cout << "can't write lagr set\n";
+
+ rval = enforce_convexity(&mb, lagrange_set);
+ if (MB_SUCCESS != rval)
+ return 1;
+
+ rval = mb.write_file("lagr.h5m", 0, 0, &lagrange_set, 1);
+ if (MB_SUCCESS != rval)
+ std::cout << "can't write lagr set\n";
+
+ Intx2MeshOnSphere worker(&mb);
+
+ double radius = 1.; // input
+
+ worker.SetRadius(radius);
+
+ //worker.SetEntityType(MBQUAD);
+
+ worker.SetErrorTolerance(gtol);
+ std::cout << "error tolerance epsilon_1=" << gtol << "\n";
+
+ EntityHandle outputSet;
+ rval = mb.create_meshset(MESHSET_SET, outputSet);
+ if (MB_SUCCESS != rval)
+ return 1;
+ rval = worker.intersect_meshes(lagrange_set, euler_set, outputSet);
+ if (MB_SUCCESS != rval)
+ return 1;
+
+ std::string opts_write("");
+ std::stringstream outf;
+ outf << "intersect1" << ".h5m";
+ rval = mb.write_file(outf.str().c_str(), 0, 0, &outputSet, 1);
+ if (MB_SUCCESS != rval)
+ std::cout << "can't write output\n";
+ double intx_area = area_on_sphere_lHuiller(&mb, outputSet, radius);
+ double arrival_area = area_on_sphere_lHuiller(&mb, euler_set, radius);
+ std::cout << " Arrival area: " << arrival_area
+ << " intersection area:" << intx_area << " rel error: "
+ << fabs((intx_area - arrival_area) / arrival_area) << "\n";
+ if (MB_SUCCESS != rval)
+ return 1;
+
+ return 0;
+}
+
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/1aadd2f328cd/
Changeset: 1aadd2f328cd
Branch: master
User: iulian07
Date: 2013-06-20 00:39:02
Summary: change the Positive and Negative tags names to red and blue parents.
Again, red is for "initial", nice, euler mesh, arrival mesh
blue is for manufactured, deformed (could be concave elements), lagrange, departure
Affected #: 4 files
diff --git a/tools/mbcslam/Intx2Mesh.cpp b/tools/mbcslam/Intx2Mesh.cpp
index a2d52e6..4c9c5d5 100644
--- a/tools/mbcslam/Intx2Mesh.cpp
+++ b/tools/mbcslam/Intx2Mesh.cpp
@@ -74,11 +74,11 @@ void Intx2Mesh::createTags()
int defaultInt = 0;
- rval = mb->tag_get_handle("Positive", 1, MB_TYPE_INTEGER, redParentTag,
+ rval = mb->tag_get_handle("RedParent", 1, MB_TYPE_INTEGER, redParentTag,
MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt);
ERRORV(rval, "can't create positive tag");
- rval = mb->tag_get_handle("Negative", 1, MB_TYPE_INTEGER, blueParentTag,
+ rval = mb->tag_get_handle("BlueParent", 1, MB_TYPE_INTEGER, blueParentTag,
MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt);
ERRORV(rval, "can't create negative tag");
@@ -156,17 +156,17 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
createTags(); //
EntityHandle startBlue, startRed;
- Range rs1;
- Range rs2;
mb->get_entities_by_dimension(mbs1, 2, rs1);
mb->get_entities_by_dimension(mbs2, 2, rs2);
- while (!rs2.empty())
+ Range rs22=rs2; // a copy of the initial range; we will remove from it elements as we
+ // advance ; rs2 is needed for marking the polygon to the red parent
+ while (!rs22.empty())
{
for (Range::iterator it = rs1.begin(); it != rs1.end(); it++)
{
startBlue = *it;
int found = 0;
- for (Range::iterator it2 = rs2.begin(); it2 != rs2.end() && !found; it2++)
+ for (Range::iterator it2 = rs22.begin(); it2 != rs22.end() && !found; it2++)
{
startRed = *it2;
double area = 0;
@@ -260,12 +260,11 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
int nb[MAXEDGES] = {0};
int nr[MAXEDGES] = {0};
- int nsidesBlue; // nsides red is the same, nsidesBlue changes for next element
- // area is in 2d, points are in 3d (on a sphere), back-projected
- // intersection points (could include the vertices of initial quads)
- // means no intersection on the side (markers)
- // means no intersection on the side (markers)
- // nc [j] = 1 means that the side j (from j to j+1) of blue poly intersects the
+ int nsidesBlue; ///
+ // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
+ // intersection points could include the vertices of initial elements
+ // nb [j] = 0 means no intersection on the side k for element blue (markers)
+ // nb [j] = 1 means that the side j (from j to j+1) of blue poly intersects the
// red poly. A potential next poly is the red poly that is adjacent to this side
computeIntersectionBetweenRedAndBlue(/* red */currentRed, blueT, P, nP,
area, nb, nr, nsidesBlue, nsidesRed);
@@ -339,8 +338,8 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
}
} // end while (!localBlue.empty())
- // here, we are finished with redCurrent, take it out of the rs2 range (red, arrival mesh)
- rs2.erase(currentRed);
+ // here, we are finished with redCurrent, take it out of the rs22 range (red, arrival mesh)
+ rs22.erase(currentRed);
// also, look at its neighbors, and add to the seeds a next one
EntityHandle redNeighbors[MAXEDGES];
diff --git a/tools/mbcslam/Intx2Mesh.hpp b/tools/mbcslam/Intx2Mesh.hpp
index b0424d9..8f5947d 100644
--- a/tools/mbcslam/Intx2Mesh.hpp
+++ b/tools/mbcslam/Intx2Mesh.hpp
@@ -95,6 +95,8 @@ protected: // so it can be accessed in derived classes, InPlane and OnSphere
EntityHandle mbs1;
EntityHandle mbs2;
+ Range rs1;// range set 1 (departure set, lagrange set, blue set, manufactured set)
+ Range rs2;// range set 2 (arrival set, euler set, red set, initial set)
EntityHandle outSet; // will contain intersection
diff --git a/tools/mbcslam/Intx2MeshInPlane.cpp b/tools/mbcslam/Intx2MeshInPlane.cpp
index 4978f91..496903a 100644
--- a/tools/mbcslam/Intx2MeshInPlane.cpp
+++ b/tools/mbcslam/Intx2MeshInPlane.cpp
@@ -340,10 +340,10 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, int nsRed, EntityHandle blue,
mb->create_element(MBPOLYGON, foundIds, nP, polyNew);
mb->add_entities(outSet, &polyNew, 1);
- // tag it with the ids from red and blue
- int id = mb->id_from_handle(blue);
+ // tag it with the index ids from red and blue sets
+ int id = rs1.index(blue); // index starts from 0
mb->tag_set_data(blueParentTag, &polyNew, 1, &id);
- id = mb->id_from_handle(red);
+ id = rs2.index(red);
mb->tag_set_data(redParentTag, &polyNew, 1, &id);
static int count=0;
diff --git a/tools/mbcslam/Intx2MeshOnSphere.cpp b/tools/mbcslam/Intx2MeshOnSphere.cpp
index c8590ad..4fea220 100644
--- a/tools/mbcslam/Intx2MeshOnSphere.cpp
+++ b/tools/mbcslam/Intx2MeshOnSphere.cpp
@@ -342,10 +342,10 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, int nsRed, EntityHandle blue,
mb->create_element(MBPOLYGON, foundIds, nP, polyNew);
mb->add_entities(outSet, &polyNew, 1);
- // tag it with the ids from red and blue
- int id = mb->id_from_handle(blue);
+ // tag it with the index ids from red and blue sets
+ int id = rs1.index(blue); // index starts from 0
mb->tag_set_data(blueParentTag, &polyNew, 1, &id);
- id = mb->id_from_handle(red);
+ id = rs2.index(red);
mb->tag_set_data(redParentTag, &polyNew, 1, &id);
static int count=0;
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list