[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Jan 6 13:51:14 CST 2014
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/1491a201d7f2/
Changeset: 1491a201d7f2
Branch: None
User: iulian07
Date: 2014-01-06 19:12:08
Summary: robustness issues with intersection
appeared at diffusion example at time step 5;
2 very close vertices, on departure and arrival meshes
were left outside everything:
intersection of edges is very strict (interval [0, 1]
in a normalized segment)
and interior points in a convex poly is also too strict
(comparing areas with 0.)
this fix relaxes the interior issues in convex polygon
by introducing an epsilon square (epsilon_area) which
is set at the same time with distance epsilon (epsilon_1)
so area check is now area < -epsilon_area
still need to think about interior segment normalization
Affected #: 7 files
diff --git a/tools/mbcslam/CslamUtils.cpp b/tools/mbcslam/CslamUtils.cpp
index 6cca02d..e806d61 100644
--- a/tools/mbcslam/CslamUtils.cpp
+++ b/tools/mbcslam/CslamUtils.cpp
@@ -30,7 +30,7 @@ double area2D(double *a, double *b, double *c)
// (b-a)x(c-a) / 2
return ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2;
}
-int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES])
+int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES], double epsilon_area)
{
// 2 triangles, 3 corners, is the corner of X in Y?
// Y must have a positive area
@@ -54,7 +54,7 @@ int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int
double area2 = (B[0] - A[0]) * (C[1] - A[1])
- (C[0] - A[0]) * (B[1] - A[1]);
- if (area2 < 0.)
+ if (area2 < -epsilon_area)
{
inside = 0;
break;
diff --git a/tools/mbcslam/CslamUtils.hpp b/tools/mbcslam/CslamUtils.hpp
index 05775df..66f4e99 100644
--- a/tools/mbcslam/CslamUtils.hpp
+++ b/tools/mbcslam/CslamUtils.hpp
@@ -21,7 +21,7 @@ namespace moab
{
double dist2(double * a, double * b);
double area2D(double *a, double *b, double *c);
-int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES]);
+int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES], double epsilon_area);
int SortAndRemoveDoubles2(double * P, int & nP, double epsilon);
// the marks will show what edges of blue intersect the red
diff --git a/tools/mbcslam/Intx2Mesh.cpp b/tools/mbcslam/Intx2Mesh.cpp
index c1da3fd..625a442 100644
--- a/tools/mbcslam/Intx2Mesh.cpp
+++ b/tools/mbcslam/Intx2Mesh.cpp
@@ -18,7 +18,7 @@
namespace moab {
-Intx2Mesh::Intx2Mesh(Interface * mbimpl):mb(mbimpl), parcomm(NULL), myTree(NULL), remote_cells(NULL)
+Intx2Mesh::Intx2Mesh(Interface * mbimpl): mb(mbimpl), parcomm(NULL), remote_cells(NULL)
{
dbg_1=0;
box_error=0;
diff --git a/tools/mbcslam/Intx2Mesh.hpp b/tools/mbcslam/Intx2Mesh.hpp
index 2daa052..d54b0b5 100644
--- a/tools/mbcslam/Intx2Mesh.hpp
+++ b/tools/mbcslam/Intx2Mesh.hpp
@@ -34,7 +34,6 @@ namespace moab {
// forward declarations
class ParallelComm;
-class AdaptiveKDTree;
class TupleList;
class Intx2Mesh
@@ -68,7 +67,7 @@ public:
ErrorCode GetOrderedNeighbors(EntityHandle set, EntityHandle quad,
EntityHandle neighbors[MAXEDGES]);
- void SetErrorTolerance(double eps) { epsilon_1=eps;}
+ void SetErrorTolerance(double eps) { epsilon_1=eps; epsilon_area = eps*eps/2;}
//void SetEntityType (EntityType tp) { type=tp;}
@@ -97,7 +96,8 @@ public:
ErrorCode correct_intersection_points_positions();
- void enable_debug() {dbg_1=1;};
+ void enable_debug() {dbg_1 = 1;}
+ void disable_debug() {dbg_1 = 0;}
protected: // so it can be accessed in derived classes, InPlane and OnSphere
Interface * mb;
@@ -139,10 +139,10 @@ protected: // so it can be accessed in derived classes, InPlane and OnSphere
std::vector<std::vector<EntityHandle> *> extraNodesVec;
double epsilon_1;
+ double epsilon_area;
ParallelComm * parcomm;
- AdaptiveKDTree *myTree;
std::vector<double> allBoxes;
double box_error;
/* \brief Local root of the kdtree */
diff --git a/tools/mbcslam/Intx2MeshInPlane.cpp b/tools/mbcslam/Intx2MeshInPlane.cpp
index 9fb2a94..2ddea3a 100644
--- a/tools/mbcslam/Intx2MeshInPlane.cpp
+++ b/tools/mbcslam/Intx2MeshInPlane.cpp
@@ -106,7 +106,7 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
}
int side[MAXEDGES] = { 0 };// this refers to what side? blue or red?
- int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side);
+ int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsBlue; k++)
@@ -133,7 +133,7 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
}
nP += extraPoints;
- extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side);
+ extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsRed; k++)
diff --git a/tools/mbcslam/Intx2MeshOnSphere.cpp b/tools/mbcslam/Intx2MeshOnSphere.cpp
index a1d786b..ef4d219 100644
--- a/tools/mbcslam/Intx2MeshOnSphere.cpp
+++ b/tools/mbcslam/Intx2MeshOnSphere.cpp
@@ -124,7 +124,7 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
return 1; // some unforeseen error
int side[MAXEDGES] = { 0 };// this refers to what side? blue or red?
- int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side);
+ int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsBlue; k++)
@@ -144,7 +144,7 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
}
nP += extraPoints;
- extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side);
+ extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side, epsilon_area);
if (extraPoints >= 1)
{
for (int k = 0; k < nsRed; k++)
@@ -370,10 +370,11 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, int nsRed, EntityHandle blue,
std::stringstream fff;
fff << "file0" << count<< ".vtk";
- mb->write_mesh(fff.str().c_str(), &outSet, 1);
+ mb->write_mesh(fff.str().c_str(), &outSet, 1);
}
}
+ disable_debug();
delete[] foundIds;
foundIds = NULL;
return 0;
diff --git a/tools/mbcslam/intx_on_sphere_test.cpp b/tools/mbcslam/intx_on_sphere_test.cpp
index 350a39b..360aa28 100644
--- a/tools/mbcslam/intx_on_sphere_test.cpp
+++ b/tools/mbcslam/intx_on_sphere_test.cpp
@@ -69,10 +69,12 @@ int main(int argc, char* argv[])
return 1;
Intx2MeshOnSphere worker(mb);
- double radius= 6. * sqrt(3.) / 2; // input
- worker.SetErrorTolerance(radius*1.e-8);
+
+ worker.SetErrorTolerance(R*1.e-8);
//worker.SetEntityType(moab::MBQUAD);
- worker.SetRadius(radius);
+ worker.SetRadius(R);
+ //worker.enable_debug();
+
rval = worker.intersect_meshes(sf1, sf2, outputSet);
//compute total area with 2 methods
@@ -81,11 +83,15 @@ int main(int argc, char* argv[])
rval = mb->write_mesh(newFile, &outputSet, 1);
if (MB_SUCCESS != rval)
return 1;
- double area_method1 = area_on_sphere_lHuiller(mb, outputSet, radius);
- double area_method2 = area_on_sphere(mb, outputSet, radius);
+ double initial_area = area_on_sphere_lHuiller(mb, sf1, R);
+ double area_method1 = area_on_sphere_lHuiller(mb, outputSet, R);
+ double area_method2 = area_on_sphere(mb, outputSet, R);
+ std::cout << "initial area: " << initial_area << "\n";
std::cout<< " area with l'Huiller: " << area_method1 << " with Girard: " << area_method2<< "\n";
- std::cout << " relative difference: " << fabs(area_method1-area_method2)/area_method1 << "\n";
+ std::cout << " relative difference areas " << fabs(area_method1-area_method2)/area_method1 << "\n";
+ std::cout << " relative error " << fabs(area_method1-initial_area)/area_method1 << "\n";
+
return 0;
}
https://bitbucket.org/fathomteam/moab/commits/ec2239392ca2/
Changeset: ec2239392ca2
Branch: master
User: iulian07
Date: 2014-01-06 20:44:27
Summary: modify imesh example
create the euler mesh before calling the wrapper
so the wrapper takes as input the euler mesh set
that contains the arrival mesh; departure point
is defined as DP tag
A tracer field is also a tag (TracerAverage), a value for
each cell
After calling the wrapper from C (intx_imesh) or from Fortran
(advection), the tracer is updated and saved to a new file)
All this works in parallel too
with
mpiexec -np 2 intx_imesh
or
mpiexec -np 2 advection
(need to prepare the input file using create_dp )
Affected #: 3 files
diff --git a/tools/mbcslam/advection.F90 b/tools/mbcslam/advection.F90
index fb102ba..ad92f88 100644
--- a/tools/mbcslam/advection.F90
+++ b/tools/mbcslam/advection.F90
@@ -26,7 +26,7 @@ program advection
use ISO_C_BINDING
implicit none
iMesh_Instance, INTENT(IN) , VALUE :: instance
- iBase_EntityHandle, INTENT(OUT) :: opEulerSet
+ iBase_EntityHandle, INTENT(IN), VALUE :: opEulerSet
integer(c_int) , INTENT (OUT) :: ierr
END SUBROUTINE update_tracer
END INTERFACE
@@ -35,18 +35,22 @@ program advection
! imesh is the instance handle
iMesh_Instance imesh
- ! ents, verts will be arrays storing vertex/entity handles
- iBase_EntityHandle, pointer :: ents, verts
+ ! cells will be storing 2d cells
+ TYPE(C_PTR) cells_ptr
+ iBase_EntityHandle, pointer :: cells(:)
+ INTEGER ents_alloc, ents_size
+
iBase_EntitySetHandle root_set
iBase_EntitySetHandle opEulerSet
CHARACTER (LEN=200) options
CHARACTER (LEN=200) filename
CHARACTER (LEN=200) optionswrite
CHARACTER (LEN=200) outname
- TYPE(C_PTR) :: vertsPtr, entsPtr
+ ! TYPE(C_PTR) :: vertsPtr, entsPtr
integer rank, sz, ierr
integer lenopt, lenname
+ integer isList
#ifdef USE_MPI
! local variables for parallel runs
@@ -95,6 +99,24 @@ program advection
%VAL(lenopt) ) !119) !strlen(options));
CHECK("fail to load mesh in parallel ")
+ isList = 0
+ call iMesh_createEntSet(%VAL(imesh), %VAL(isList), opEulerSet, ierr)
+ CHECK("fail to create euler set ")
+
+ ents_alloc = 0
+ ents_size = 0
+
+ call iMesh_getEntities(%VAL(imesh),%VAL(root_set), &
+ %VAL(iBase_FACE), &
+ %VAL(iMesh_ALL_TOPOLOGIES), cells_ptr, &
+ ents_alloc, ents_size, ierr);
+ CHECK("fail to get 2d entities ")
+
+ call c_f_pointer(cells_ptr, cells, [ents_size])
+
+ call iMesh_addEntArrToSet(%VAL(imesh), cells, %VAL(ents_size), &
+ %VAL(opEulerSet), ierr)
+
call update_tracer(imesh, opEulerSet, ierr)
CHECK("fail to update tracer ")
diff --git a/tools/mbcslam/intx_imesh.cpp b/tools/mbcslam/intx_imesh.cpp
index 0c49ac4..cb7b494 100644
--- a/tools/mbcslam/intx_imesh.cpp
+++ b/tools/mbcslam/intx_imesh.cpp
@@ -16,7 +16,7 @@
#ifdef __cplusplus
extern "C" {
#endif
- void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr);
+ void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle eulerSet, int * ierr);
#ifdef __cplusplus
} // extern "C"
@@ -71,7 +71,23 @@ int main(int argc, char* argv[]){
printf("There's %d entity sets here on process rank %d \n", num_sets, rank);
iBase_EntitySetHandle euler_set;
- update_tracer( imesh, &euler_set, &ierr);
+
+ iMesh_createEntSet(imesh, 0, &euler_set, &ierr);
+ IMESH_ASSERT(ierr);
+
+ iBase_EntityHandle *cells = NULL;
+ int ents_alloc = 0;
+ int ents_size = 0;
+
+ iMesh_getEntities(imesh, root, iBase_FACE, iMesh_ALL_TOPOLOGIES, &cells,
+ &ents_alloc, &ents_size, &ierr);
+ IMESH_ASSERT(ierr);
+
+ iMesh_addEntArrToSet(imesh, cells, ents_size, euler_set, &ierr);
+
+ IMESH_ASSERT(ierr);
+
+ update_tracer( imesh, euler_set, &ierr);
IMESH_ASSERT(ierr);
// write everything
diff --git a/tools/mbcslam/wrap_intx.cpp b/tools/mbcslam/wrap_intx.cpp
index e38c897..b973acf 100644
--- a/tools/mbcslam/wrap_intx.cpp
+++ b/tools/mbcslam/wrap_intx.cpp
@@ -23,25 +23,14 @@ bool debug = false;
extern "C" {
#endif
-void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr)
+void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, int * ierr)
{
Range ents;
moab::Interface * mb =MOABI;
*ierr =1;
- ErrorCode rval = mb->get_entities_by_dimension(0, 2, ents);// root 0
- ERRORV(rval, "can't get all 2d entities from root");
- EntityHandle euler_set;
- //EntityHandle lagr_set;
-
- rval = mb->create_meshset(MESHSET_SET, euler_set);
- ERRORV(rval , "can't create arrival mesh set");
-
- *opEulerSet = (iBase_EntitySetHandle)euler_set;
-
- rval = mb->add_entities(euler_set, ents);
- ERRORV(rval , "can't add ents to arrival set");
+ EntityHandle euler_set = (EntityHandle) imesh_euler_set;
Intx2MeshOnSphere worker(mb);
worker.SetRadius(radius);
@@ -49,7 +38,8 @@ void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet,
worker.SetErrorTolerance(gtol);
EntityHandle covering_lagr_set;
- rval = mb->create_meshset(MESHSET_SET, covering_lagr_set);
+
+ ErrorCode rval = mb->create_meshset(MESHSET_SET, covering_lagr_set);
ERRORV(rval , "can't create covering set ");
// we need to update the correlation tag and remote tuples
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