From 06f2b772d525db91bc7b71d69ffb1acca210355f Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 26 Jul 2018 13:36:04 -0400 Subject: [PATCH 01/19] phasta: std for cerr --- phasta/phiotimer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phasta/phiotimer.cc b/phasta/phiotimer.cc index e55dcde04..69f4e88dd 100644 --- a/phasta/phiotimer.cc +++ b/phasta/phiotimer.cc @@ -44,7 +44,7 @@ static size_t phastaio_getCyclesPerMicroSec() { cycles = t1 - t0; cpus = ((double)cycles)/(usec); if(!PCU_Comm_Self()) - cerr << "cycles " << cycles << " us " << usec + std::cerr << "cycles " << cycles << " us " << usec << " cycles per micro second " << cpus << "\n"; return cpus; } From 8629587735f1eb8b3185a1e7f4a7d4f866cc88b4 Mon Sep 17 00:00:00 2001 From: jared Date: Wed, 25 Jul 2018 10:57:16 -0400 Subject: [PATCH 02/19] introduce reduceFieldData and ReductionOp implements #159 --- apf/apf.cc | 9 ++++++-- apf/apf.h | 55 ++++++++++++++++++++++++++++++++++++++++++++- apf/apfFieldData.cc | 16 ++++++------- apf/apfFieldData.h | 4 +++- pumi/pumi_field.cc | 2 +- 5 files changed, 72 insertions(+), 14 deletions(-) diff --git a/apf/apf.cc b/apf/apf.cc index 13eca372a..263ab7bf0 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -424,11 +424,16 @@ void synchronize(Field* f, Sharing* shr) synchronizeFieldData(f->getData(), shr); } -void accumulate(Field* f, Sharing* shr) +void accumulate(Field* f, Sharing* shr, bool delete_shr) { - accumulateFieldData(f->getData(), shr); + sharedReduction(f, shr, delete_shr, ReductionSum() ); } +void sharedReduction(Field* f, Sharing* shr, bool delete_shr, + const ReductionOp& reduce_op ) +{ + reduceFieldData(f->getData(), shr, delete_shr, reduce_op); +} void fail(const char* why) { fprintf(stderr,"APF FAILED: %s\n",why); diff --git a/apf/apf.h b/apf/apf.h index 5921afeaf..a1b7c2bd9 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -37,6 +37,50 @@ class VectorElement; typedef VectorElement MeshElement; class FieldShape; struct Sharing; +template class ReductionOp; +template class ReductionSum; + +/** \brief Base class for applying operations to make a Field consistent + * in parallel + * \details This function gets applied pairwise to the Field values + * from every partition, resulting in a single unique value. No guarantees + * are made about the order in which this function is applied to the + * values. + */ +template +class ReductionOp +{ + public: + /* \brief apply operation, returning a single value */ + virtual T apply(T val1, T val2) const = 0; +}; + +template +class ReductionSum : public ReductionOp +{ + T apply(T val1, T val2) const { return val1 + val2; }; +}; + +template +class ReductionMin : public ReductionOp +{ + T apply(T val1, T val2) const { return ( (val1 < val2) ? val1 : val2 ); }; +}; + +template +class ReductionMax : public ReductionOp +{ + T apply(T val1, T val2) const { return ( (val1 < val2) ? val2 : val1 ); }; +}; + + +/* instantiate (is this necessary with the global consts below?) */ +template class ReductionSum; +template class ReductionMin; +template class ReductionMax; + + + /** \brief Destroys an apf::Mesh. * @@ -625,7 +669,16 @@ void synchronize(Field* f, Sharing* shr = 0); all copies of an entity and assign the sum as the value for all copies. */ -void accumulate(Field* f, Sharing* shr = 0); +void accumulate(Field* f, Sharing* shr = 0, bool delete_shr = false); + +/** \brief Apply a reduction operator along a partition boundary + \details Using the copies described by an apf::Sharing object, applies + the specified operation pairwise to the values of the field on each + partition. No guarantee is made about the order of the pairwise + application. + */ +void sharedReduction(Field* f, Sharing* shr = 0, bool delete_shr=false, + const ReductionOp& reduce_op = ReductionSum()); /** \brief Declare failure of code inside APF. \details This function prints the string as an APF diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index b89b087dc..b5705b86b 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -81,7 +81,7 @@ template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); -void accumulateFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr) +void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, const ReductionOp& reduce_op /* =ReductionSum() */) { FieldBase* f = data->getField(); Mesh* m = f->getMesh(); @@ -100,17 +100,14 @@ void accumulateFieldData(FieldDataOf* data, Sharing* shr, bool delete_sh PCU_Comm_Begin(); while ((e = m->iterate(it))) { - if (( ! data->hasEntity(e)) || m->isGhost(e) || - (shr->isOwned(e))) - continue; /* non-owners send to owners */ + if (( ! data->hasEntity(e)) || m->isGhost(e) ) + continue; /* send to all parts that can see this entity */ CopyArray copies; shr->getCopies(e, copies); int n = f->countValuesOn(e); NewArray values(n); data->get(e,&(values[0])); - /* actually, non-owners send to all others, - since apf::Sharing doesn't identify the owner */ for (size_t i = 0; i < copies.getSize(); ++i) { PCU_COMM_PACK(copies[i].peer, copies[i].entity); @@ -131,11 +128,12 @@ void accumulateFieldData(FieldDataOf* data, Sharing* shr, bool delete_sh PCU_Comm_Unpack(&(inValues[0]),n*sizeof(double)); data->get(e,&(values[0])); for (int i = 0; i < n; ++i) - values[i] += inValues[i]; + values[i] = reduce_op.apply(values[i], inValues[i]); data->set(e,&(values[0])); } - } /* broadcast back out to non-owners */ - synchronizeFieldData(data, shr, delete_shr); + } + // every partition did the reduction, so no need to broadcast result + if (delete_shr) delete shr; } template diff --git a/apf/apfFieldData.h b/apf/apfFieldData.h index 24be9cd0d..bffb0d790 100644 --- a/apf/apfFieldData.h +++ b/apf/apfFieldData.h @@ -11,6 +11,7 @@ #include #include "apfField.h" #include "apfShape.h" +#include "apf.h" // needed for ReductionOp namespace apf { @@ -35,7 +36,7 @@ class FieldDataOf; template void synchronizeFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr=false); -void accumulateFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr=false); +void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr=false, const ReductionOp& reduce_op=ReductionSum()); template void copyFieldData(FieldDataOf* from, FieldDataOf* to); @@ -57,6 +58,7 @@ class FieldDataOf : public FieldData int getElementData(MeshEntity* entity, NewArray& data); }; + } //namespace apf #endif diff --git a/pumi/pumi_field.cc b/pumi/pumi_field.cc index fa721df0d..ada73fe90 100644 --- a/pumi/pumi_field.cc +++ b/pumi/pumi_field.cc @@ -135,7 +135,7 @@ void pumi_field_synchronize(pField f, pOwnership o) void pumi_field_accumulate(pField f, pOwnership o) { - apf::accumulateFieldData(f->getData(), o, false); + apf::reduceFieldData(f->getData(), o, false); } void pumi_field_freeze(pField f) From 4ac36a5e145bdbe88d4b6f7cbeecd7bebd6fc6fc Mon Sep 17 00:00:00 2001 From: jared Date: Tue, 31 Jul 2018 13:48:25 -0400 Subject: [PATCH 03/19] test passes --- apf/apfFieldData.cc | 26 ++++++++++++++++++++++++-- apf/apfFieldData.h | 2 ++ test/CMakeLists.txt | 1 + 3 files changed, 27 insertions(+), 2 deletions(-) diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index b5705b86b..d2928145c 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -83,6 +83,12 @@ template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, const ReductionOp& reduce_op /* =ReductionSum() */) { + char fname[255]; + std::sprintf(fname, "debug2_%d.txt", PCU_Comm_Self()); + std::cout << "opening file " << fname << std::endl; + std::ofstream fout; + fout.open(fname); + fout << "Reducing field data 2" << std::endl; FieldBase* f = data->getField(); Mesh* m = f->getMesh(); FieldShape* s = f->getShape(); @@ -91,15 +97,20 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c shr = getSharing(m); delete_shr=true; } + + fout << "apf::Sharing* = " << shr << std::endl; for (int d=0; d < 4; ++d) { if ( ! s->hasNodesIn(d)) continue; + + fout << "d = " << d << std::endl; MeshEntity* e; MeshIterator* it = m->begin(d); PCU_Comm_Begin(); while ((e = m->iterate(it))) { + fout << "considering entity " << e << ", hasEntity = " << data->hasEntity(e) << ", isGhost = " << m->isGhost(e) << std::endl; if (( ! data->hasEntity(e)) || m->isGhost(e) ) continue; /* send to all parts that can see this entity */ @@ -108,8 +119,15 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c int n = f->countValuesOn(e); NewArray values(n); data->get(e,&(values[0])); + + if (copies.getSize() > 0) + fout << "\nentity = " << e << std::endl; + else + fout << "no copies, skipping" << std::endl; + for (size_t i = 0; i < copies.getSize(); ++i) { + fout << "sending entity " << copies[i].entity << " to peer " << copies[i].peer << " with value " << values[0] << std::endl; PCU_COMM_PACK(copies[i].peer, copies[i].entity); PCU_Comm_Pack(copies[i].peer, &(values[0]), n*sizeof(double)); } @@ -118,20 +136,24 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c PCU_Comm_Send(); while (PCU_Comm_Listen()) while ( ! PCU_Comm_Unpacked()) - { /* receive and add. we only care about correctness - on the owners */ + { /* receive and apply reduction */ MeshEntity* e; PCU_COMM_UNPACK(e); int n = f->countValuesOn(e); NewArray values(n); NewArray inValues(n); PCU_Comm_Unpack(&(inValues[0]),n*sizeof(double)); + fout << "receiving entity " << e << " with value " << inValues[0] << std::endl; data->get(e,&(values[0])); for (int i = 0; i < n; ++i) + { values[i] = reduce_op.apply(values[i], inValues[i]); + } data->set(e,&(values[0])); } } + + fout << "finished reduction" << std::endl; // every partition did the reduction, so no need to broadcast result if (delete_shr) delete shr; } diff --git a/apf/apfFieldData.h b/apf/apfFieldData.h index bffb0d790..9b476545b 100644 --- a/apf/apfFieldData.h +++ b/apf/apfFieldData.h @@ -13,6 +13,8 @@ #include "apfShape.h" #include "apf.h" // needed for ReductionOp +#include // DEBUGGING +#include namespace apf { class FieldData diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 946be8aba..b2148ad78 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -165,6 +165,7 @@ test_exe_func(poisson poisson.cc) test_exe_func(ph_adapt ph_adapt.cc) test_exe_func(assert_timing assert_timing.cc) test_exe_func(create_mis create_mis.cc) +test_exe_func(fieldReduce fieldReduce.cc) if(ENABLE_DSP) test_exe_func(graphdist graphdist.cc) test_exe_func(moving moving.cc) From 1e135a2f397e0b0a9a3f4963d62b9227646761e6 Mon Sep 17 00:00:00 2001 From: jared Date: Tue, 31 Jul 2018 14:08:02 -0400 Subject: [PATCH 04/19] add test to CTest --- test/testing.cmake | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/testing.cmake b/test/testing.cmake index 30765d4af..90d460e20 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -327,6 +327,11 @@ if(ENABLE_ZOLTAN) "pipe_4_.smb" "tet.smb") endif() +mpi_test(fieldReduce 4 + ./fieldReduce + "${MDIR}/pipe.${GXT}" + "pipe_4_.smb") + set(MDIR ${MESHES}/torus) mpi_test(reorder 4 ./reorder From 83fcca0a354bd0179ff0aed18a9da08ac9457b6e Mon Sep 17 00:00:00 2001 From: jared Date: Tue, 31 Jul 2018 14:11:50 -0400 Subject: [PATCH 05/19] remove debug output --- apf/apfFieldData.cc | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index d2928145c..46052920a 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -83,12 +83,6 @@ template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, const ReductionOp& reduce_op /* =ReductionSum() */) { - char fname[255]; - std::sprintf(fname, "debug2_%d.txt", PCU_Comm_Self()); - std::cout << "opening file " << fname << std::endl; - std::ofstream fout; - fout.open(fname); - fout << "Reducing field data 2" << std::endl; FieldBase* f = data->getField(); Mesh* m = f->getMesh(); FieldShape* s = f->getShape(); @@ -98,19 +92,16 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c delete_shr=true; } - fout << "apf::Sharing* = " << shr << std::endl; for (int d=0; d < 4; ++d) { if ( ! s->hasNodesIn(d)) continue; - fout << "d = " << d << std::endl; MeshEntity* e; MeshIterator* it = m->begin(d); PCU_Comm_Begin(); while ((e = m->iterate(it))) { - fout << "considering entity " << e << ", hasEntity = " << data->hasEntity(e) << ", isGhost = " << m->isGhost(e) << std::endl; if (( ! data->hasEntity(e)) || m->isGhost(e) ) continue; /* send to all parts that can see this entity */ @@ -120,14 +111,8 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c NewArray values(n); data->get(e,&(values[0])); - if (copies.getSize() > 0) - fout << "\nentity = " << e << std::endl; - else - fout << "no copies, skipping" << std::endl; - for (size_t i = 0; i < copies.getSize(); ++i) { - fout << "sending entity " << copies[i].entity << " to peer " << copies[i].peer << " with value " << values[0] << std::endl; PCU_COMM_PACK(copies[i].peer, copies[i].entity); PCU_Comm_Pack(copies[i].peer, &(values[0]), n*sizeof(double)); } @@ -143,7 +128,6 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c NewArray values(n); NewArray inValues(n); PCU_Comm_Unpack(&(inValues[0]),n*sizeof(double)); - fout << "receiving entity " << e << " with value " << inValues[0] << std::endl; data->get(e,&(values[0])); for (int i = 0; i < n; ++i) { @@ -153,7 +137,6 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c } } - fout << "finished reduction" << std::endl; // every partition did the reduction, so no need to broadcast result if (delete_shr) delete shr; } From 23d1cdc4b467b0d27cc827fd2e6c9739ea676971 Mon Sep 17 00:00:00 2001 From: jared Date: Wed, 1 Aug 2018 09:10:16 -0400 Subject: [PATCH 06/19] add missing file --- test/fieldReduce.cc | 197 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 197 insertions(+) create mode 100644 test/fieldReduce.cc diff --git a/test/fieldReduce.cc b/test/fieldReduce.cc new file mode 100644 index 000000000..03ec45d68 --- /dev/null +++ b/test/fieldReduce.cc @@ -0,0 +1,197 @@ +#include +#include +#include +#include +#include +#include +#include +#ifdef HAVE_SIMMETRIX +#include +#include +#include +#include +#endif +#include +#include +#include +#include +namespace { + +const char* modelFile = 0; +const char* meshFile = 0; +const char* outFile = 0; +int myrank = -1; +int commsize = -1; + +double getValue(apf::Vector3& coords, double addval) +{ + return coords.x() + coords.y() + 1 + addval; +} + +apf::Field* getTestField(apf::Mesh* m, const char* fname, double addval) +{ + apf::FieldShape* fshape = apf::getLagrange(1); + apf::Field* f = apf::createLagrangeField(m, fname, apf::SCALAR, 1); + + // write known data to field, check that the sum is n times the number of + // copies + apf::MeshEntity* e; + apf::MeshIterator* it; + apf::Vector3 coords; + + for (int d=0; d <= m->getDimension(); ++d) + { + if (!fshape->hasNodesIn(d)) + continue; + + it = m->begin(d); + while ( (e = m->iterate(it)) ) + { + m->getPoint(e, 0, coords); + double val = getValue(coords, addval); + apf::setScalar(f, e, 0, val); + } // end while + } // end for + + return f; +} + +bool testReduce(apf::Mesh* m, int casenum) +{ + double addval = 0; + if (casenum == 0) // sum + addval = 0; + else + addval = myrank; + + char fname[256]; + sprintf(fname, "test%d", casenum); + apf::Field* f = getTestField(m, fname, addval); + apf::FieldShape* fshape = apf::getShape(f); + apf::Sharing* shr = apf::getSharing(m); + + if (casenum == 0) + apf::accumulate(f, shr); + else if (casenum == 1) + apf::sharedReduction(f, shr, false, apf::ReductionMax()); + else if (casenum == 2) + apf::sharedReduction(f, shr, false, apf::ReductionMin()); + + // verify the result is n times the number of copies + apf::MeshEntity* e; + apf::MeshIterator* it; + apf::Vector3 coords; + bool failflag = false; + + for (int d=0; d <= m->getDimension(); ++d) + { + if (!fshape->hasNodesIn(d)) + continue; + + it = m->begin(d); + while ( (e = m->iterate(it)) ) + { + apf::CopyArray copies; // need to redeclare this every iteration so + // the size gets reset + shr->getCopies(e, copies); + int ntimes = copies.getSize() + 1; + + if (ntimes > 1) + { + // this entity has remotes, get some info + m->getPoint(e, 0, coords); + double val_f = apf::getScalar(f, e, 0); + + // get max and min peer + int maxpeer = myrank; + int minpeer = myrank; + for (std::size_t j = 0; j < copies.getSize(); ++j) + { + if (copies[j].peer > maxpeer) + maxpeer = copies[j].peer; + + if (copies[j].peer < minpeer) + minpeer = copies[j].peer; + } + + // do the test + if (casenum == 0) + { + double val = getValue(coords, addval); + failflag = ( failflag || ( std::fabs(ntimes*val - val_f) > 1e-13 ) ); + } else if (casenum == 1) // max + { + double val = getValue(coords, maxpeer); + failflag = ( failflag || ( std::fabs(val - val_f) > 1e-13 ) ); + } else if (casenum == 2) { // min + double val = getValue(coords, minpeer); + failflag = ( failflag || ( std::fabs(val - val_f) > 1e-13 ) ); + } // end if casenum + + } // end if ntimes + } // end while + } // end for + + return failflag; +} + +void freeMesh(apf::Mesh* m) +{ + m->destroyNative(); + apf::destroyMesh(m); +} + +void getConfig(int argc, char** argv) +{ + if ( argc != 3 ) { + if ( !PCU_Comm_Self() ) + printf("Usage: %s \n", argv[0]); + MPI_Finalize(); + exit(EXIT_FAILURE); + } + modelFile = argv[1]; + meshFile = argv[2]; + outFile = argv[3]; + PCU_ALWAYS_ASSERT(commsize == 4); +} + +} + +int main(int argc, char** argv) +{ + MPI_Init(&argc,&argv); + PCU_Comm_Init(); + MPI_Comm_rank(PCU_Get_Comm(), &myrank); + MPI_Comm_size(PCU_Get_Comm(), &commsize); +#ifdef HAVE_SIMMETRIX + MS_init(); + SimModel_start(); + Sim_readLicenseFile(0); + gmi_sim_start(); + gmi_register_sim(); +#endif + gmi_register_mesh(); + getConfig(argc,argv); + gmi_model* g = 0; + g = gmi_load(modelFile); + apf::Mesh2* m = 0; + m = apf::loadMdsMesh(g, meshFile); + + bool failflag = false; + for (int i=0; i < 3; ++i) + failflag = failflag || testReduce(m, i); + + + freeMesh(m); +#ifdef HAVE_SIMMETRIX + gmi_sim_stop(); + Sim_unregisterAllKeys(); + SimModel_stop(); + MS_exit(); +#endif + PCU_Comm_Free(); + MPI_Finalize(); + + return failflag; +} + From 082d1ad084d00272a0bfec344c771a336734238e Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 1 Aug 2018 10:11:35 -0400 Subject: [PATCH 07/19] generate supports surface or volume only mesh generation --- test/generate.cc | 144 +++++++++++++++++++++++++++++++++------------ test/testing.cmake | 14 ++++- 2 files changed, 119 insertions(+), 39 deletions(-) diff --git a/test/generate.cc b/test/generate.cc index ecf2c0be5..d36a3724c 100644 --- a/test/generate.cc +++ b/test/generate.cc @@ -17,6 +17,9 @@ #include #include +#include //cout +#include //option parser + #ifdef SIM_PARASOLID #include "SimParasolidKrnl.h" #endif @@ -37,8 +40,12 @@ pAManager SModel_attManager(pModel model); namespace { +int should_log = 0; +int disable_volume = 0; +int disable_surface = 0; std::string modelFile; std::string nativeModelFile; +std::string surfaceMeshFile; std::string caseName; std::string outMeshFile; @@ -76,56 +83,117 @@ pParMesh generate(pGModel mdl, std::string meshCaseName) { MS_processSimModelerAdvMeshingAtts(mcaseFile, mcase); AttCase_setModel(mcase, mdl); - pParMesh pmesh = PM_new(0, mdl, PMU_size()); + pParMesh pmesh; + if( ! surfaceMeshFile.empty() && + disable_surface && !disable_volume ) { + //load the surface mesh instead of creating it + pmesh = PM_load(surfaceMeshFile.c_str(), mdl, NULL); + PM_setTotalNumParts(pmesh, PMU_size()); //enable parallel volume meshing + } else { + //create an empty surface mesh + pmesh = PM_new(0, mdl, PMU_size()); + } - const double stime = MPI_Wtime(); - if(0==PCU_Comm_Self()) { - printf("Meshing surface..."); fflush(stdout); + if( !disable_surface ) { + const double stime = MPI_Wtime(); + if(0==PCU_Comm_Self()) { + printf("Meshing surface..."); fflush(stdout); + } + pSurfaceMesher surfaceMesher = SurfaceMesher_new(mcase, pmesh); + SurfaceMesher_execute(surfaceMesher, NULL); + SurfaceMesher_delete(surfaceMesher); + if(0==PCU_Comm_Self()) + printf(" %f seconds\n", MPI_Wtime()-stime); + if( ! surfaceMeshFile.empty() ) { + if(0==PCU_Comm_Self()) + printf(" writing surface mesh %s\n", surfaceMeshFile.c_str()); + PM_write(pmesh, surfaceMeshFile.c_str(), NULL); + } } - pSurfaceMesher surfaceMesher = SurfaceMesher_new(mcase, pmesh); - SurfaceMesher_execute(surfaceMesher, NULL); - SurfaceMesher_delete(surfaceMesher); - if(0==PCU_Comm_Self()) - printf(" %f seconds\n", MPI_Wtime()-stime); - const double vtime = MPI_Wtime(); - if(0==PCU_Comm_Self()) { - printf("Meshing volume..."); fflush(stdout); + if( !disable_volume ) { + const double vtime = MPI_Wtime(); + if(0==PCU_Comm_Self()) { + printf("Meshing volume..."); fflush(stdout); + } + pVolumeMesher volumeMesher = VolumeMesher_new(mcase, pmesh); + VolumeMesher_execute(volumeMesher, NULL); + VolumeMesher_delete(volumeMesher); + if(0==PCU_Comm_Self()) + printf(" %f seconds\n", MPI_Wtime()-vtime); } - pVolumeMesher volumeMesher = VolumeMesher_new(mcase, pmesh); - VolumeMesher_execute(volumeMesher, NULL); - VolumeMesher_delete(volumeMesher); - if(0==PCU_Comm_Self()) - printf(" %f seconds\n", MPI_Wtime()-vtime); return pmesh; } -void getConfig(int argc, char** argv) -{ - if (argc < 3) { - if(0==PCU_Comm_Self()) { - printf("Usage: %s ", argv[0]); - printf(" to generate a mesh on a GeomSim model\n"); - printf(" or: %s \n", argv[0]); - printf(" to generate a mesh using the specified case name using the SimModeler" - "model which references the native parasolid or acis model\n"); +void getConfig(int argc, char** argv) { + opterr = 0; + + static struct option long_opts[] = { + {"enable-log", no_argument, &should_log, 1}, + {"disable-volume", no_argument, &disable_volume, 1}, + {"disable-surface", no_argument, &disable_surface, 1}, + {"native-model", required_argument, 0, 'n'}, + {"surface-mesh", required_argument, 0, 'm'}, + {0, 0, 0, 0} // terminate the option array + }; + + const char* usage="" + "[options] \n" + "options:\n" + " --enable-log Enable Simmetrix logging\n" + " --disable-volume Disable volume mesh generation\n" + " --disable-surface Disable suface mesh generation\n" + " --native-model=/path/to/model Load the native Parasolid or ACIS model that the GeomSim model uses\n"; + + nativeModelFile = ""; + surfaceMeshFile = ""; + int option_index = 0; + while(1) { + int c = getopt_long(argc, argv, "", long_opts, &option_index); + if (c == -1) break; //end of options + switch (c) { + case 0: // enable-log|disable-volume|disable-surf + if (!PCU_Comm_Self()) + printf ("read arg %d\n", c); + break; + case 'n': + nativeModelFile = std::string(optarg); + break; + case 'm': + surfaceMeshFile = std::string(optarg); + break; + case '?': + if (!PCU_Comm_Self()) + printf ("warning: skipping unrecognized option\n"); + break; + default: + if (!PCU_Comm_Self()) + printf("Usage %s %s", argv[0], usage); + exit(EXIT_FAILURE); } - MPI_Finalize(); - exit(EXIT_SUCCESS); } - modelFile = argv[1]; - if (argc == 3) { - nativeModelFile = ""; - outMeshFile = caseName = argv[2]; - } else if (argc == 4) { - nativeModelFile = argv[2]; - outMeshFile = caseName = argv[3]; + + + if(argc-optind != 2) { + if (!PCU_Comm_Self()) + printf("Usage %s %s", argv[0], usage); + exit(EXIT_FAILURE); } + int i=optind; + modelFile = std::string(argv[i++]); + outMeshFile = caseName = std::string(argv[i++]); outMeshFile.append("/"); - if(0==PCU_Comm_Self()) { - printf("Inputs: model \'%s\' native model \'%s\' case \'%s\'\n", - modelFile.c_str(), nativeModelFile.c_str(), caseName.c_str()); + + if (!PCU_Comm_Self()) { + std::cout << "enable_log " << should_log << + " disable_surface " << disable_surface << + " disable_volume " << disable_volume << + " native-model " << nativeModelFile << + " model " << modelFile << + " surface mesh " << surfaceMeshFile << + " case name " << caseName << + " output mesh" << outMeshFile << "\n"; } } diff --git a/test/testing.cmake b/test/testing.cmake index 90d460e20..8df14c230 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -511,11 +511,23 @@ if(ENABLE_SIMMETRIX) ./generate "${MDIR}/upright.smd" "67k") + mpi_test(parallel_meshgen_surf 4 + ./generate + "--disable-volume" + "--surface-mesh=${MDIR}/67k_surf.sms" + "${MDIR}/upright.smd" + "67k") + mpi_test(parallel_meshgen_vol 4 + ./generate + "--disable-surface" + "--surface-mesh=${MDIR}/67k_surf_ref.sms" + "${MDIR}/upright.smd" + "67k") if(SIM_PARASOLID) mpi_test(parallel_meshgen_para 4 ./generate + "--native-model=${MDIR}/upright.x_t" "${MDIR}/upright.smd" - "${MDIR}/upright.x_t" "67k") endif() # adapt_meshgen uses the output of parallel_meshgen From 686a9232ef1db43fea94c65df41953c65861e011 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Tue, 31 Jul 2018 16:47:33 -0400 Subject: [PATCH 08/19] implement FieldData::clone() this adds an implementation to clone field data which is needed for apf::convert to copy fields. closes #162 --- apf/apf.cc | 7 ++++ apf/apf.h | 18 +++++--- apf/apfArrayData.cc | 8 ++++ apf/apfCoordData.cc | 8 ++++ apf/apfCoordData.h | 1 + apf/apfFieldData.cc | 5 --- apf/apfFieldData.h | 3 +- apf/apfUserData.cc | 9 ++++ apf/apfUserData.h | 8 +++- test/CMakeLists.txt | 1 + test/testing.cmake | 1 + test/verify_convert.cc | 95 ++++++++++++++++++++++++++++++++++++++++++ 12 files changed, 151 insertions(+), 13 deletions(-) create mode 100644 test/verify_convert.cc diff --git a/apf/apf.cc b/apf/apf.cc index 263ab7bf0..ead73b867 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -468,6 +468,13 @@ Field* createUserField(Mesh* m, const char* name, int valueType, FieldShape* s, return makeField(m, name, valueType, 0, s, new UserData(f)); } +void updateUserField(Field* field, Function* newFunc) +{ + UserData* ud = dynamic_cast(field->getData()); + // ud will be null if the data is not user data + if (ud) ud->setFunction(newFunc); +} + void copyData(Field* to, Field* from) { copyFieldData(from->getData(), to->getData()); diff --git a/apf/apf.h b/apf/apf.h index a1b7c2bd9..4e8153133 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -721,14 +721,22 @@ struct Function virtual void eval(MeshEntity* e, double* result) = 0; }; -/** \brief Create a Field from a user's analytic function. - \details This field will use no memory, has values on all - nodes, and calls the user Function for all value queries. - Writing to this field does nothing. - */ +/* \brief Create a Field from a user's analytic function. + * \details This field will use no memory, has values on all + * nodes, and calls the user Function for all value queries. + * Writing to this field does nothing. + * \warning if you copy the mesh with apf::convert you may want to use + * apf::updateUserField to update function that this field refers to. This is + * extremely important if the analytic function you use references user fields. + */ Field* createUserField(Mesh* m, const char* name, int valueType, FieldShape* s, Function* f); +/* \brief update the analytic function from a user field + * \details this field updates the apf::Function which the UserField refers to + */ +void updateUserField(Field* field, Function* newFunc); + /** \brief Compute a nodal gradient field from a nodal input field \details given a nodal field, compute approximate nodal gradient values by giving each node a volume-weighted average of the diff --git a/apf/apfArrayData.cc b/apf/apfArrayData.cc index 289afaa28..3d76f7c0a 100644 --- a/apf/apfArrayData.cc +++ b/apf/apfArrayData.cc @@ -79,6 +79,14 @@ class ArrayDataOf : public FieldDataOf T* getDataArray() { return this->dataArray; } + virtual FieldData* clone() { + //FieldData* newData = new TagDataOf(); + FieldData* newData = new ArrayDataOf(); + newData->init(this->field); + copyFieldData(static_cast*>(newData), + static_cast*>(this->field->getData())); + return newData; + } private: /* data variables go here */ diff --git a/apf/apfCoordData.cc b/apf/apfCoordData.cc index 0578540e5..b74031c82 100644 --- a/apf/apfCoordData.cc +++ b/apf/apfCoordData.cc @@ -11,6 +11,14 @@ namespace apf { +FieldData* CoordData::clone() { + FieldData* newData = new CoordData(); + newData->init(field); + copyFieldData(static_cast*>(newData), + static_cast*>(field->getData())); + return newData; + +} void CoordData::init(FieldBase* f) { FieldData::field = f; diff --git a/apf/apfCoordData.h b/apf/apfCoordData.h index 002428f9b..8866b79f1 100644 --- a/apf/apfCoordData.h +++ b/apf/apfCoordData.h @@ -33,6 +33,7 @@ class CoordData : public FieldDataOf virtual void get(MeshEntity* e, double* data); virtual void set(MeshEntity* e, double const* data); virtual bool isFrozen() { return false; } + virtual FieldData* clone(); private: Mesh* mesh; FieldShape* shape; diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index 46052920a..a6133bc84 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -10,11 +10,6 @@ FieldData::~FieldData() { } -FieldData* FieldData::clone() -{ - abort(); -} - void FieldData::rename(const char*) { abort(); diff --git a/apf/apfFieldData.h b/apf/apfFieldData.h index 9b476545b..8d28014b6 100644 --- a/apf/apfFieldData.h +++ b/apf/apfFieldData.h @@ -25,7 +25,7 @@ class FieldData virtual bool hasEntity(MeshEntity* e) = 0; virtual void removeEntity(MeshEntity* e) = 0; virtual bool isFrozen() = 0; - virtual FieldData* clone(); + virtual FieldData* clone() = 0; virtual void rename(const char* newName); FieldBase* getField() {return field;} protected: @@ -58,6 +58,7 @@ class FieldDataOf : public FieldData void setNodeComponents(MeshEntity* e, int node, T const* components); void getNodeComponents(MeshEntity* e, int node, T* components); int getElementData(MeshEntity* entity, NewArray& data); + virtual FieldData* clone()=0; }; diff --git a/apf/apfUserData.cc b/apf/apfUserData.cc index e324f189c..15ee6a310 100644 --- a/apf/apfUserData.cc +++ b/apf/apfUserData.cc @@ -42,4 +42,13 @@ bool UserData::isFrozen() return false; } +FieldData* UserData::clone() +{ + FieldData* newData = new UserData(function); + newData->init(field); + copyFieldData(static_cast*>(newData), + static_cast*>(field->getData())); + return newData; +} + } diff --git a/apf/apfUserData.h b/apf/apfUserData.h index 4802cd822..96980727e 100644 --- a/apf/apfUserData.h +++ b/apf/apfUserData.h @@ -21,11 +21,15 @@ struct UserData : public FieldDataOf void get(MeshEntity* e, double* data); void set(MeshEntity* e, double const* data); bool isFrozen(); + virtual FieldData* clone(); + // using const * const gives an error on gcc/7.3.0 because the return is an + // r-value which cannot be modified anyways + Function const* getFunction() const { return function; } + void setFunction(Function* func) { function = func; } + private: Function* function; }; } #endif - - diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b2148ad78..d06168999 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,6 +20,7 @@ util_exe_func(describe describe.cc) test_exe_func(quality quality.cc) test_exe_func(writeVtxPtn writeVtxPtn.cc) test_exe_func(verify_2nd_order_shapes verify_2nd_order_shapes.cc) +test_exe_func(verify_convert verify_convert.cc) # Geometric model utilities if(ENABLE_SIMMETRIX) diff --git a/test/testing.cmake b/test/testing.cmake index 8df14c230..582ea920b 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -33,6 +33,7 @@ mpi_test(integrate 1 ./integrate) mpi_test(qr_test 1 ./qr) mpi_test(base64 1 ./base64) mpi_test(tensor_test 1 ./tensor) +mpi_test(verify_convert 1 ./verify_convert) if(ENABLE_SIMMETRIX) mpi_test(in_closure_of 1 diff --git a/test/verify_convert.cc b/test/verify_convert.cc new file mode 100644 index 000000000..72302603a --- /dev/null +++ b/test/verify_convert.cc @@ -0,0 +1,95 @@ +/* + * this test verifies that the convert process properly clones the underlying + * fields + */ +#include +#include +#include +#include +#include +#include +#include +#include +apf::Mesh2* createEmptyMesh() +{ + gmi_model* mdl = gmi_load(".null"); + return apf::makeEmptyMdsMesh(mdl, 3, false); +} +apf::Mesh2* createMesh() +{ + apf::Mesh2* m = createEmptyMesh(); + apf::MeshEntity* verts[2]; + verts[0] = + m->createVertex(NULL, apf::Vector3(0, 0, 0), apf::Vector3(0, 0, 0)); + verts[1] = + m->createVertex(NULL, apf::Vector3(1, 0, 0), apf::Vector3(1, 0, 0)); + m->createEntity(apf::Mesh::Type::EDGE, NULL, verts); + return m; +} +class twox : public apf::Function { + private: + apf::Field* x; + + public: + virtual void eval(apf::MeshEntity* e, double* result) + { + result[0] = 2 * apf::getScalar(x, e, 0); + } + twox(apf::Field* x) : x(x) {} +}; +int main(int argc, char* argv[]) +{ + MPI_Init(&argc, &argv); + PCU_Comm_Init(); + gmi_register_null(); + apf::Mesh* m1 = createMesh(); + apf::Mesh2* m2 = createEmptyMesh(); + // create field on m1 + apf::Field* f = apf::createLagrangeField(m1, "field1", apf::SCALAR, 1); + apf::Function* func = new twox(f); + apf::Field* uf = + apf::createUserField(m1, "ufield1", apf::SCALAR, apf::getShape(f), func); + // loop over all vertices in mesh and set values + int count = 1; + apf::MeshIterator* it = m1->begin(0); + while (apf::MeshEntity* vert = m1->iterate(it)) { + apf::setScalar(f, vert, 0, count++); + } + m1->end(it); + // verify the user field works on mesh 1 + it = m1->begin(0); + while (apf::MeshEntity* vert = m1->iterate(it)) { + double val = apf::getScalar(f, vert, 0); + double uval = apf::getScalar(uf, vert, 0); + if (!(std::abs(uval - 2 * double(val)) < 1E-15)) return 1; + } + m1->end(it); + // copy m1 to m2 + apf::convert(m1, m2); + apf::Field* f2 = m2->findField("field1"); + apf::Field* uf2 = m2->findField("ufield1"); + // update the user field to reference the field in mesh 2 + apf::updateUserField(uf2, new twox(f2)); + count = 1; + it = m2->begin(0); + while (apf::MeshEntity* vert = m2->iterate(it)) { + double val = apf::getScalar(f2, vert, 0); + double uval = apf::getScalar(uf2, vert, 0); + assert(std::abs(val - count) < 1E-15); + assert(std::abs(uval - 2 * double(count)) < 1E-15); + apf::setScalar(f2, vert, 0, 18.0); + // make sure that the function updated properly + uval = apf::getScalar(uf2, vert, 0); + assert(std::abs(uval - 36.0) < 1E-15); + ++count; + } + m2->end(it); + delete func; + m1->destroyNative(); + m2->destroyNative(); + apf::destroyMesh(m1); + apf::destroyMesh(m2); + PCU_Comm_Free(); + MPI_Finalize(); + return 0; +} From 9d7d3537ec4e48caaebe3b738bcedebaa8685fe2 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Thu, 2 Aug 2018 00:54:54 -0400 Subject: [PATCH 09/19] update field of numbering when converted/cloned --- apf/apfConvert.cc | 28 ++++++++++++++++++++++++---- apf/apfNumbering.cc | 10 +++++++++- apf/apfNumbering.h | 8 ++++++++ test/verify_convert.cc | 38 ++++++++++++++++++++++++++++++++++++-- 4 files changed, 77 insertions(+), 7 deletions(-) diff --git a/apf/apfConvert.cc b/apf/apfConvert.cc index fc4b163cc..bea0370b5 100644 --- a/apf/apfConvert.cc +++ b/apf/apfConvert.cc @@ -222,8 +222,18 @@ class Converter { for (int i = 0; i < inMesh->countNumberings(); ++i) { Numbering* in = inMesh->getNumbering(i); - Numbering* out = createNumbering(outMesh, - getName(in), getShape(in), countComponents(in)); + Numbering* out; + if (getField(in)) { + // here we assume that the fields have already been copied into the + // mesh + Field* outField = outMesh->findField(getName(getField(in))); + PCU_DEBUG_ASSERT(outField); + out = createNumbering(outField); + } + else { + out = createNumbering(outMesh, getName(in), getShape(in), + countComponents(in)); + } convertNumbering(in, out); } } @@ -231,8 +241,18 @@ class Converter { for (int i = 0; i < inMesh->countGlobalNumberings(); ++i) { GlobalNumbering* in = inMesh->getGlobalNumbering(i); - GlobalNumbering* out = createGlobalNumbering(outMesh, - getName(in), getShape(in), countComponents(in)); + GlobalNumbering* out; + if (getField(in)) { + // here we assume that the fields have already been copied into the + // mesh + Field* outField = outMesh->findField(getName(getField(in))); + PCU_DEBUG_ASSERT(outField); + out = createGlobalNumbering(outField); + } + else { + out = createGlobalNumbering(outMesh, getName(in), getShape(in), + countComponents(in)); + } convertGlobalNumbering(in, out); } } diff --git a/apf/apfNumbering.cc b/apf/apfNumbering.cc index 1a70b2272..8b05d6c3a 100644 --- a/apf/apfNumbering.cc +++ b/apf/apfNumbering.cc @@ -461,6 +461,14 @@ GlobalNumbering* createGlobalNumbering( return n; } +GlobalNumbering* createGlobalNumbering(Field* f) +{ + GlobalNumbering* n = new GlobalNumbering(); + n->init(f); + f->getMesh()->addGlobalNumbering(n); + return n; +} + FieldShape* getShape(GlobalNumbering* n) { return n->getShape(); @@ -601,5 +609,5 @@ void getNodes(GlobalNumbering* n, DynamicArray& nodes) getFieldNodes(n,nodes); } +Field* getField(GlobalNumbering* n) { return n->getField(); } } - diff --git a/apf/apfNumbering.h b/apf/apfNumbering.h index 35b02c334..f4f339e4d 100644 --- a/apf/apfNumbering.h +++ b/apf/apfNumbering.h @@ -159,6 +159,11 @@ GlobalNumbering* createGlobalNumbering( const char* name, FieldShape* shape, int components=1); + +/** \brief Create a Numbering of degrees of freedom of a Field. + */ +GlobalNumbering* createGlobalNumbering(Field* f); + FieldShape* getShape(GlobalNumbering* n); const char* getName(GlobalNumbering* n); /** \brief get the mesh associated with a global numbering */ @@ -176,6 +181,9 @@ long getNumber(GlobalNumbering* n, MeshEntity* e, int node, int component=0); /** \brief get an element's global node numbers */ int getElementNumbers(GlobalNumbering* n, MeshEntity* e, NewArray& numbers); +/** \brief get the field being numbered + */ +Field* getField(GlobalNumbering* n); /** \brief converts a local numbering into a global numbering. \param destroy Should the input Numbering* be destroyed? diff --git a/test/verify_convert.cc b/test/verify_convert.cc index 72302603a..8569b431b 100644 --- a/test/verify_convert.cc +++ b/test/verify_convert.cc @@ -1,12 +1,13 @@ /* * this test verifies that the convert process properly clones the underlying - * fields + * fields and numberings */ #include #include #include #include #include +#include #include #include #include @@ -49,11 +50,25 @@ int main(int argc, char* argv[]) apf::Function* func = new twox(f); apf::Field* uf = apf::createUserField(m1, "ufield1", apf::SCALAR, apf::getShape(f), func); + // create numbering and global numbering with and without fields to make sure + // they transfer properly + apf::Numbering* numWithField = apf::createNumbering(f); + apf::Numbering* numNoField = apf::createNumbering( + m1, "noField", apf::getShape(f), apf::countComponents(f)); + apf::GlobalNumbering* globalNumWithField = apf::createGlobalNumbering(f); + apf::GlobalNumbering* globalNumNoField = apf::createGlobalNumbering( + m1, "noField", apf::getShape(f), apf::countComponents(f)); // loop over all vertices in mesh and set values int count = 1; apf::MeshIterator* it = m1->begin(0); while (apf::MeshEntity* vert = m1->iterate(it)) { - apf::setScalar(f, vert, 0, count++); + apf::setScalar(f, vert, 0, count); + // set the numberings + apf::number(numWithField, vert, 0, 0, count); + apf::number(numNoField, vert, 0, 0, count); + apf::number(globalNumWithField, vert, 0, count); + apf::number(globalNumNoField, vert, 0, count); + ++count; } m1->end(it); // verify the user field works on mesh 1 @@ -68,6 +83,20 @@ int main(int argc, char* argv[]) apf::convert(m1, m2); apf::Field* f2 = m2->findField("field1"); apf::Field* uf2 = m2->findField("ufield1"); + + apf::Numbering* numWithField2 = m2->findNumbering(apf::getName(numWithField)); + apf::Numbering* numNoField2 = m2->findNumbering(apf::getName(numNoField)); + apf::GlobalNumbering* globalNumWithField2 = + m2->findGlobalNumbering(apf::getName(globalNumWithField)); + apf::GlobalNumbering* globalNumNoField2 = + m2->findGlobalNumbering(apf::getName(globalNumNoField)); + // all of these numberings should exist + assert(numWithField2 && numNoField2 && globalNumWithField2 && + globalNumNoField2); + // make sure the fields of the numberings match up properly as they should be + // the new field copied into the new mesh + assert(getField(numWithField2) == f2); + assert(getField(globalNumWithField2) == f2); // update the user field to reference the field in mesh 2 apf::updateUserField(uf2, new twox(f2)); count = 1; @@ -81,6 +110,11 @@ int main(int argc, char* argv[]) // make sure that the function updated properly uval = apf::getScalar(uf2, vert, 0); assert(std::abs(uval - 36.0) < 1E-15); + // check to make sure the numberings have the correct values + assert(getNumber(numWithField2, vert, 0, 0) == count); + assert(getNumber(numNoField2, vert, 0, 0) == count); + assert(getNumber(globalNumWithField2, vert, 0, 0) == count); + assert(getNumber(globalNumNoField2, vert, 0, 0) == count); ++count; } m2->end(it); From 2c216475946b0b1e18f2f7a4df3877c300da618a Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Thu, 2 Aug 2018 14:37:15 -0400 Subject: [PATCH 10/19] clone tag data in convert --- apf/apfConvert.cc | 81 ++++++++++++++++++++++++++++++++++++++++++ test/verify_convert.cc | 15 +++++++- 2 files changed, 95 insertions(+), 1 deletion(-) diff --git a/apf/apfConvert.cc b/apf/apfConvert.cc index bea0370b5..6f30f2863 100644 --- a/apf/apfConvert.cc +++ b/apf/apfConvert.cc @@ -31,6 +31,9 @@ class Converter convertFields(); convertNumberings(); convertGlobalNumberings(); + // this must be called after anything that might create tags e.g. fields + // or numberings to avoid problems with tag duplication + convertTags(); outMesh->acceptChanges(); } ModelEntity* getNewModelFromOld(ModelEntity* oldC) @@ -210,6 +213,48 @@ class Converter } } } + void convertTag(Mesh* inMesh, MeshTag* in, Mesh* outMesh, MeshTag* out) + { + for (int d = 0; d <= 3; ++d) { + int tagType = inMesh->getTagType(in); + int tagSize = inMesh->getTagSize(in); + PCU_DEBUG_ASSERT(tagType == outMesh->getTagType(out)); + PCU_DEBUG_ASSERT(tagSize == outMesh->getTagSize(out)); + MeshIterator* it = inMesh->begin(d); + MeshEntity* e; + while ((e = inMesh->iterate(it))) { + if(inMesh->hasTag(e, in)) { + // these initializations cannot go into the cases due to compiler + // warnings on gcc 7.3.0 + double* dblData; + int* intData; + long* lngData; + switch (tagType) { + case apf::Mesh::TagType::DOUBLE: + dblData = new double[tagSize]; + inMesh->getDoubleTag(e, in, dblData); + outMesh->setDoubleTag(newFromOld[e], out, dblData); + break; + case apf::Mesh::TagType::INT: + intData = new int[tagSize]; + inMesh->getIntTag(e, in, intData); + outMesh->setIntTag(newFromOld[e], out, intData); + break; + case apf::Mesh::TagType::LONG: + lngData = new long[tagSize]; + inMesh->getLongTag(e, in, lngData); + outMesh->setLongTag(newFromOld[e], out, lngData); + break; + default: + std::cerr << "Tried to convert unknown tag type\n"; + abort(); + break; + } + } + } + inMesh->end(it); + } + } void convertFields() { for (int i = 0; i < inMesh->countFields(); ++i) { @@ -256,6 +301,42 @@ class Converter convertGlobalNumbering(in, out); } } + void convertTags() + { + DynamicArray tags; + inMesh->getTags(tags); + for (std::size_t i = 0; i < tags.getSize(); ++i) { + apf::MeshTag* in = tags[i]; + PCU_DEBUG_ASSERT(in); + // create a new tag on the outMesh + int tagType = inMesh->getTagType(in); + int tagSize = inMesh->getTagSize(in); + const char* tagName = inMesh->getTagName(in); + PCU_DEBUG_ASSERT(tagName); + // need to make sure that the tag wasn't already created by a field or + // numbering + if (!outMesh->findTag(tagName)) { + apf::MeshTag* out = NULL; + switch (tagType) { + case apf::Mesh::TagType::DOUBLE: + out = outMesh->createDoubleTag(tagName, tagSize); + break; + case apf::Mesh::TagType::INT: + out = outMesh->createIntTag(tagName, tagSize); + break; + case apf::Mesh::TagType::LONG: + out = outMesh->createLongTag(tagName, tagSize); + break; + default: + std::cerr << "Tried to convert unknown tag type\n"; + abort(); + } + PCU_DEBUG_ASSERT(out); + // copy the tag on the inMesh to the outMesh + convertTag(inMesh, in, outMesh, out); + } + } + } void convertQuadratic() { if (inMesh->getShape() != getLagrange(2) && inMesh->getShape() != getSerendipity()) diff --git a/test/verify_convert.cc b/test/verify_convert.cc index 8569b431b..73d7e86e1 100644 --- a/test/verify_convert.cc +++ b/test/verify_convert.cc @@ -1,6 +1,6 @@ /* * this test verifies that the convert process properly clones the underlying - * fields and numberings + * fields, numberings, and tags */ #include #include @@ -58,6 +58,8 @@ int main(int argc, char* argv[]) apf::GlobalNumbering* globalNumWithField = apf::createGlobalNumbering(f); apf::GlobalNumbering* globalNumNoField = apf::createGlobalNumbering( m1, "noField", apf::getShape(f), apf::countComponents(f)); + // create an integer tag + apf::MeshTag* intTag = m1->createIntTag("intTag", 1); // loop over all vertices in mesh and set values int count = 1; apf::MeshIterator* it = m1->begin(0); @@ -68,6 +70,9 @@ int main(int argc, char* argv[]) apf::number(numNoField, vert, 0, 0, count); apf::number(globalNumWithField, vert, 0, count); apf::number(globalNumNoField, vert, 0, count); + // set the tag + // int tagData[1] = {count}; + m1->setIntTag(vert, intTag, &count); ++count; } m1->end(it); @@ -99,6 +104,9 @@ int main(int argc, char* argv[]) assert(getField(globalNumWithField2) == f2); // update the user field to reference the field in mesh 2 apf::updateUserField(uf2, new twox(f2)); + // find the copied tag data + apf::MeshTag* intTag2 = m2->findTag("intTag"); + assert(intTag2); count = 1; it = m2->begin(0); while (apf::MeshEntity* vert = m2->iterate(it)) { @@ -115,6 +123,11 @@ int main(int argc, char* argv[]) assert(getNumber(numNoField2, vert, 0, 0) == count); assert(getNumber(globalNumWithField2, vert, 0, 0) == count); assert(getNumber(globalNumNoField2, vert, 0, 0) == count); + // check that the correct tag data was recovered + int* data = new int[1]; + m2->getIntTag(vert, intTag2, data); + assert(*data == count); + delete[] data; ++count; } m2->end(it); From 6eee54768b0a0929ea063f53e2300d16fb23a6d3 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 2 Aug 2018 22:06:06 -0400 Subject: [PATCH 11/19] generate: move logging, use flag --- test/generate.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/generate.cc b/test/generate.cc index d36a3724c..760133f21 100644 --- a/test/generate.cc +++ b/test/generate.cc @@ -261,6 +261,8 @@ pNativeModel loadNativeModel() { void simStart() { SimModel_start(); SimPartitionedMesh_start(NULL,NULL); + if(should_log) + Sim_logOn("generate_sim.log"); MS_init(); SimModel_start(); #ifdef SIM_PARASOLID @@ -300,7 +302,6 @@ int main(int argc, char** argv) getConfig(argc,argv); simStart(); - Sim_logOn("generate_sim.log"); pNativeModel nm = loadNativeModel(); pGModel simModel = GM_load(modelFile.c_str(), nm, NULL); From 597243d18454e8a31317c0350a1d8512c322211f Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 3 Aug 2018 14:42:04 -0400 Subject: [PATCH 12/19] generate: add surface mesh to usage stmt --- test/generate.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/generate.cc b/test/generate.cc index 760133f21..09f78824d 100644 --- a/test/generate.cc +++ b/test/generate.cc @@ -141,10 +141,11 @@ void getConfig(int argc, char** argv) { const char* usage="" "[options] \n" "options:\n" - " --enable-log Enable Simmetrix logging\n" - " --disable-volume Disable volume mesh generation\n" - " --disable-surface Disable suface mesh generation\n" - " --native-model=/path/to/model Load the native Parasolid or ACIS model that the GeomSim model uses\n"; + " --enable-log Enable Simmetrix logging\n" + " --disable-volume Disable volume mesh generation\n" + " --disable-surface Disable suface mesh generation\n" + " --native-model=/path/to/model Load the native Parasolid or ACIS model that the GeomSim model uses\n" + " --surfaceMeshFile=/path/to/surfaceMesh read or write the surface mesh - depends on generation mode\n"; nativeModelFile = ""; surfaceMeshFile = ""; From 938d7f97c957e6b1dfc4770ced96f6a28ac6a71a Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 3 Aug 2018 14:46:09 -0400 Subject: [PATCH 13/19] generate: fix usage stmt --- test/generate.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/generate.cc b/test/generate.cc index 09f78824d..c9b94b93a 100644 --- a/test/generate.cc +++ b/test/generate.cc @@ -145,7 +145,7 @@ void getConfig(int argc, char** argv) { " --disable-volume Disable volume mesh generation\n" " --disable-surface Disable suface mesh generation\n" " --native-model=/path/to/model Load the native Parasolid or ACIS model that the GeomSim model uses\n" - " --surfaceMeshFile=/path/to/surfaceMesh read or write the surface mesh - depends on generation mode\n"; + " --surface-mesh=/path/to/surfaceMesh read or write the surface mesh - depends on generation mode\n"; nativeModelFile = ""; surfaceMeshFile = ""; From d783b19e0b879e7c0f8b18a987c1eac313f90343 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Wed, 8 Aug 2018 07:56:11 -0400 Subject: [PATCH 14/19] Revert "Merge pull request #161 from JaredCrean2/jc_generalize_accumulate2" This reverts commit 70d9f1cd794002ed3bcef710db0dd2b233f8887d, reversing changes made to 06f2b772d525db91bc7b71d69ffb1acca210355f. --- apf/apf.cc | 9 +- apf/apf.h | 55 +------------ apf/apfFieldData.cc | 25 +++--- apf/apfFieldData.h | 6 +- pumi/pumi_field.cc | 2 +- test/CMakeLists.txt | 1 - test/fieldReduce.cc | 197 -------------------------------------------- test/testing.cmake | 5 -- 8 files changed, 16 insertions(+), 284 deletions(-) delete mode 100644 test/fieldReduce.cc diff --git a/apf/apf.cc b/apf/apf.cc index 263ab7bf0..13eca372a 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -424,16 +424,11 @@ void synchronize(Field* f, Sharing* shr) synchronizeFieldData(f->getData(), shr); } -void accumulate(Field* f, Sharing* shr, bool delete_shr) +void accumulate(Field* f, Sharing* shr) { - sharedReduction(f, shr, delete_shr, ReductionSum() ); + accumulateFieldData(f->getData(), shr); } -void sharedReduction(Field* f, Sharing* shr, bool delete_shr, - const ReductionOp& reduce_op ) -{ - reduceFieldData(f->getData(), shr, delete_shr, reduce_op); -} void fail(const char* why) { fprintf(stderr,"APF FAILED: %s\n",why); diff --git a/apf/apf.h b/apf/apf.h index a1b7c2bd9..5921afeaf 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -37,50 +37,6 @@ class VectorElement; typedef VectorElement MeshElement; class FieldShape; struct Sharing; -template class ReductionOp; -template class ReductionSum; - -/** \brief Base class for applying operations to make a Field consistent - * in parallel - * \details This function gets applied pairwise to the Field values - * from every partition, resulting in a single unique value. No guarantees - * are made about the order in which this function is applied to the - * values. - */ -template -class ReductionOp -{ - public: - /* \brief apply operation, returning a single value */ - virtual T apply(T val1, T val2) const = 0; -}; - -template -class ReductionSum : public ReductionOp -{ - T apply(T val1, T val2) const { return val1 + val2; }; -}; - -template -class ReductionMin : public ReductionOp -{ - T apply(T val1, T val2) const { return ( (val1 < val2) ? val1 : val2 ); }; -}; - -template -class ReductionMax : public ReductionOp -{ - T apply(T val1, T val2) const { return ( (val1 < val2) ? val2 : val1 ); }; -}; - - -/* instantiate (is this necessary with the global consts below?) */ -template class ReductionSum; -template class ReductionMin; -template class ReductionMax; - - - /** \brief Destroys an apf::Mesh. * @@ -669,16 +625,7 @@ void synchronize(Field* f, Sharing* shr = 0); all copies of an entity and assign the sum as the value for all copies. */ -void accumulate(Field* f, Sharing* shr = 0, bool delete_shr = false); - -/** \brief Apply a reduction operator along a partition boundary - \details Using the copies described by an apf::Sharing object, applies - the specified operation pairwise to the values of the field on each - partition. No guarantee is made about the order of the pairwise - application. - */ -void sharedReduction(Field* f, Sharing* shr = 0, bool delete_shr=false, - const ReductionOp& reduce_op = ReductionSum()); +void accumulate(Field* f, Sharing* shr = 0); /** \brief Declare failure of code inside APF. \details This function prints the string as an APF diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index 46052920a..b89b087dc 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -81,7 +81,7 @@ template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); -void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, const ReductionOp& reduce_op /* =ReductionSum() */) +void accumulateFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr) { FieldBase* f = data->getField(); Mesh* m = f->getMesh(); @@ -91,26 +91,26 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c shr = getSharing(m); delete_shr=true; } - for (int d=0; d < 4; ++d) { if ( ! s->hasNodesIn(d)) continue; - MeshEntity* e; MeshIterator* it = m->begin(d); PCU_Comm_Begin(); while ((e = m->iterate(it))) { - if (( ! data->hasEntity(e)) || m->isGhost(e) ) - continue; /* send to all parts that can see this entity */ + if (( ! data->hasEntity(e)) || m->isGhost(e) || + (shr->isOwned(e))) + continue; /* non-owners send to owners */ CopyArray copies; shr->getCopies(e, copies); int n = f->countValuesOn(e); NewArray values(n); data->get(e,&(values[0])); - + /* actually, non-owners send to all others, + since apf::Sharing doesn't identify the owner */ for (size_t i = 0; i < copies.getSize(); ++i) { PCU_COMM_PACK(copies[i].peer, copies[i].entity); @@ -121,7 +121,8 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c PCU_Comm_Send(); while (PCU_Comm_Listen()) while ( ! PCU_Comm_Unpacked()) - { /* receive and apply reduction */ + { /* receive and add. we only care about correctness + on the owners */ MeshEntity* e; PCU_COMM_UNPACK(e); int n = f->countValuesOn(e); @@ -130,15 +131,11 @@ void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr, c PCU_Comm_Unpack(&(inValues[0]),n*sizeof(double)); data->get(e,&(values[0])); for (int i = 0; i < n; ++i) - { - values[i] = reduce_op.apply(values[i], inValues[i]); - } + values[i] += inValues[i]; data->set(e,&(values[0])); } - } - - // every partition did the reduction, so no need to broadcast result - if (delete_shr) delete shr; + } /* broadcast back out to non-owners */ + synchronizeFieldData(data, shr, delete_shr); } template diff --git a/apf/apfFieldData.h b/apf/apfFieldData.h index 9b476545b..24be9cd0d 100644 --- a/apf/apfFieldData.h +++ b/apf/apfFieldData.h @@ -11,10 +11,7 @@ #include #include "apfField.h" #include "apfShape.h" -#include "apf.h" // needed for ReductionOp -#include // DEBUGGING -#include namespace apf { class FieldData @@ -38,7 +35,7 @@ class FieldDataOf; template void synchronizeFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr=false); -void reduceFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr=false, const ReductionOp& reduce_op=ReductionSum()); +void accumulateFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr=false); template void copyFieldData(FieldDataOf* from, FieldDataOf* to); @@ -60,7 +57,6 @@ class FieldDataOf : public FieldData int getElementData(MeshEntity* entity, NewArray& data); }; - } //namespace apf #endif diff --git a/pumi/pumi_field.cc b/pumi/pumi_field.cc index ada73fe90..fa721df0d 100644 --- a/pumi/pumi_field.cc +++ b/pumi/pumi_field.cc @@ -135,7 +135,7 @@ void pumi_field_synchronize(pField f, pOwnership o) void pumi_field_accumulate(pField f, pOwnership o) { - apf::reduceFieldData(f->getData(), o, false); + apf::accumulateFieldData(f->getData(), o, false); } void pumi_field_freeze(pField f) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b2148ad78..946be8aba 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -165,7 +165,6 @@ test_exe_func(poisson poisson.cc) test_exe_func(ph_adapt ph_adapt.cc) test_exe_func(assert_timing assert_timing.cc) test_exe_func(create_mis create_mis.cc) -test_exe_func(fieldReduce fieldReduce.cc) if(ENABLE_DSP) test_exe_func(graphdist graphdist.cc) test_exe_func(moving moving.cc) diff --git a/test/fieldReduce.cc b/test/fieldReduce.cc deleted file mode 100644 index 03ec45d68..000000000 --- a/test/fieldReduce.cc +++ /dev/null @@ -1,197 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#ifdef HAVE_SIMMETRIX -#include -#include -#include -#include -#endif -#include -#include -#include -#include -namespace { - -const char* modelFile = 0; -const char* meshFile = 0; -const char* outFile = 0; -int myrank = -1; -int commsize = -1; - -double getValue(apf::Vector3& coords, double addval) -{ - return coords.x() + coords.y() + 1 + addval; -} - -apf::Field* getTestField(apf::Mesh* m, const char* fname, double addval) -{ - apf::FieldShape* fshape = apf::getLagrange(1); - apf::Field* f = apf::createLagrangeField(m, fname, apf::SCALAR, 1); - - // write known data to field, check that the sum is n times the number of - // copies - apf::MeshEntity* e; - apf::MeshIterator* it; - apf::Vector3 coords; - - for (int d=0; d <= m->getDimension(); ++d) - { - if (!fshape->hasNodesIn(d)) - continue; - - it = m->begin(d); - while ( (e = m->iterate(it)) ) - { - m->getPoint(e, 0, coords); - double val = getValue(coords, addval); - apf::setScalar(f, e, 0, val); - } // end while - } // end for - - return f; -} - -bool testReduce(apf::Mesh* m, int casenum) -{ - double addval = 0; - if (casenum == 0) // sum - addval = 0; - else - addval = myrank; - - char fname[256]; - sprintf(fname, "test%d", casenum); - apf::Field* f = getTestField(m, fname, addval); - apf::FieldShape* fshape = apf::getShape(f); - apf::Sharing* shr = apf::getSharing(m); - - if (casenum == 0) - apf::accumulate(f, shr); - else if (casenum == 1) - apf::sharedReduction(f, shr, false, apf::ReductionMax()); - else if (casenum == 2) - apf::sharedReduction(f, shr, false, apf::ReductionMin()); - - // verify the result is n times the number of copies - apf::MeshEntity* e; - apf::MeshIterator* it; - apf::Vector3 coords; - bool failflag = false; - - for (int d=0; d <= m->getDimension(); ++d) - { - if (!fshape->hasNodesIn(d)) - continue; - - it = m->begin(d); - while ( (e = m->iterate(it)) ) - { - apf::CopyArray copies; // need to redeclare this every iteration so - // the size gets reset - shr->getCopies(e, copies); - int ntimes = copies.getSize() + 1; - - if (ntimes > 1) - { - // this entity has remotes, get some info - m->getPoint(e, 0, coords); - double val_f = apf::getScalar(f, e, 0); - - // get max and min peer - int maxpeer = myrank; - int minpeer = myrank; - for (std::size_t j = 0; j < copies.getSize(); ++j) - { - if (copies[j].peer > maxpeer) - maxpeer = copies[j].peer; - - if (copies[j].peer < minpeer) - minpeer = copies[j].peer; - } - - // do the test - if (casenum == 0) - { - double val = getValue(coords, addval); - failflag = ( failflag || ( std::fabs(ntimes*val - val_f) > 1e-13 ) ); - } else if (casenum == 1) // max - { - double val = getValue(coords, maxpeer); - failflag = ( failflag || ( std::fabs(val - val_f) > 1e-13 ) ); - } else if (casenum == 2) { // min - double val = getValue(coords, minpeer); - failflag = ( failflag || ( std::fabs(val - val_f) > 1e-13 ) ); - } // end if casenum - - } // end if ntimes - } // end while - } // end for - - return failflag; -} - -void freeMesh(apf::Mesh* m) -{ - m->destroyNative(); - apf::destroyMesh(m); -} - -void getConfig(int argc, char** argv) -{ - if ( argc != 3 ) { - if ( !PCU_Comm_Self() ) - printf("Usage: %s \n", argv[0]); - MPI_Finalize(); - exit(EXIT_FAILURE); - } - modelFile = argv[1]; - meshFile = argv[2]; - outFile = argv[3]; - PCU_ALWAYS_ASSERT(commsize == 4); -} - -} - -int main(int argc, char** argv) -{ - MPI_Init(&argc,&argv); - PCU_Comm_Init(); - MPI_Comm_rank(PCU_Get_Comm(), &myrank); - MPI_Comm_size(PCU_Get_Comm(), &commsize); -#ifdef HAVE_SIMMETRIX - MS_init(); - SimModel_start(); - Sim_readLicenseFile(0); - gmi_sim_start(); - gmi_register_sim(); -#endif - gmi_register_mesh(); - getConfig(argc,argv); - gmi_model* g = 0; - g = gmi_load(modelFile); - apf::Mesh2* m = 0; - m = apf::loadMdsMesh(g, meshFile); - - bool failflag = false; - for (int i=0; i < 3; ++i) - failflag = failflag || testReduce(m, i); - - - freeMesh(m); -#ifdef HAVE_SIMMETRIX - gmi_sim_stop(); - Sim_unregisterAllKeys(); - SimModel_stop(); - MS_exit(); -#endif - PCU_Comm_Free(); - MPI_Finalize(); - - return failflag; -} - diff --git a/test/testing.cmake b/test/testing.cmake index 8df14c230..396cf5593 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -327,11 +327,6 @@ if(ENABLE_ZOLTAN) "pipe_4_.smb" "tet.smb") endif() -mpi_test(fieldReduce 4 - ./fieldReduce - "${MDIR}/pipe.${GXT}" - "pipe_4_.smb") - set(MDIR ${MESHES}/torus) mpi_test(reorder 4 ./reorder From 86b4cddf101af2e56477f43c6b2c170433f78ad7 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Wed, 8 Aug 2018 11:13:32 -0400 Subject: [PATCH 15/19] address compile issue brought up by @cwsmith --- apf/apfConvert.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/apf/apfConvert.cc b/apf/apfConvert.cc index 6f30f2863..ca617e06e 100644 --- a/apf/apfConvert.cc +++ b/apf/apfConvert.cc @@ -7,6 +7,7 @@ #include "apfNumbering.h" #include #include +#include namespace apf { From e7b17b609db66b15f4a3775d2957c11ddf326cc8 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Fri, 10 Aug 2018 16:43:34 -0400 Subject: [PATCH 16/19] Implemented canGetClosestPoint and isOnModel --- apf/apfMesh.cc | 15 +++++++++++++++ apf/apfMesh.h | 4 ++++ gmi/gmi.c | 5 +++++ gmi/gmi.h | 2 ++ 4 files changed, 26 insertions(+) diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc index 10c61a713..1b0c3f088 100644 --- a/apf/apfMesh.cc +++ b/apf/apfMesh.cc @@ -175,6 +175,11 @@ bool Mesh::canSnap() return gmi_can_eval(getModel()); } +bool Mesh::canGetClosestPoint() +{ + return gmi_can_get_closest_point(getModel()); +} + bool Mesh::canGetModelNormal() { return gmi_has_normal(getModel()); @@ -252,6 +257,16 @@ bool Mesh::isInClosureOf(ModelEntity* g, ModelEntity* target){ return (res == 1) ? true : false; } +bool Mesh::isOnModel(ModelEntity* g, Vector3 p, double scale) +{ + Vector3 to; + double param[2]; + gmi_ent* c = (gmi_ent*)g; + gmi_closest_point(getModel(), c, &p[0], &to[0], param); + double ratio = (to - p).getLength() / scale; + return ratio < 0.001; +} + void Mesh::getPoint(MeshEntity* e, int node, Vector3& p) { getVector(coordinateField,e,node,p); diff --git a/apf/apfMesh.h b/apf/apfMesh.h index c000f77dd..29c3631b4 100644 --- a/apf/apfMesh.h +++ b/apf/apfMesh.h @@ -285,6 +285,8 @@ class Mesh ModelEntity* findModelEntity(int type, int tag); /** \brief return true if the geometric model supports snapping */ bool canSnap(); + /** \brief return true if the geometric model supports get closest point */ + bool canGetClosestPoint(); /** \brief return true if the geometric model supports normal computation */ bool canGetModelNormal(); /** \brief evaluate parametric coordinate (p) as a spatial point (x) */ @@ -310,6 +312,8 @@ class Mesh Vector3 const& param, Vector3& x); /** \brief checks if g is in the closure of the target */ bool isInClosureOf(ModelEntity* g, ModelEntity* target); + /** \brief checks if p is on model g */ + bool isOnModel(ModelEntity* g, Vector3 p, double scale); /** \brief get the distribution of the mesh's coordinate field */ FieldShape* getShape() const; /** \brief get the mesh's coordinate field */ diff --git a/gmi/gmi.c b/gmi/gmi.c index 8ea4d6002..71f9b9efd 100644 --- a/gmi/gmi.c +++ b/gmi/gmi.c @@ -85,6 +85,11 @@ int gmi_can_eval(struct gmi_model* m) return m->ops->eval != NULL; } +int gmi_can_get_closest_point(struct gmi_model* m) +{ + return m->ops->closest_point != NULL; +} + int gmi_has_normal(struct gmi_model* m) { return m->ops->normal != NULL; diff --git a/gmi/gmi.h b/gmi/gmi.h index 548f788ad..48a4447dd 100644 --- a/gmi/gmi.h +++ b/gmi/gmi.h @@ -145,6 +145,8 @@ struct gmi_ent* gmi_find(struct gmi_model* m, int dim, int tag); struct gmi_set* gmi_adjacent(struct gmi_model* m, struct gmi_ent* e, int dim); /** \brief check whether the model implements gmi_eval */ int gmi_can_eval(struct gmi_model* m); +/** \brief check whether the model gmi_closest_point */ +int gmi_can_get_closest_point(struct gmi_model* m); /** \brief check whether the model implements gmi_normal */ int gmi_has_normal(struct gmi_model* m); /** \brief evaluate the parametric definition of a model boundary entity From 9c3069bbd1bb7c1d015be8b7e39e322c52bf01cd Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Fri, 10 Aug 2018 16:44:31 -0400 Subject: [PATCH 17/19] Quadratic simx meshes can be imported to Bezier The inversion preserves the validity of the curved elements --- crv/crvBezier.h | 2 +- crv/crvCurveMesh.cc | 21 ++++++++++++++++++--- crv/crvShapeHandler.cc | 8 ++++---- 3 files changed, 23 insertions(+), 8 deletions(-) diff --git a/crv/crvBezier.h b/crv/crvBezier.h index 97c8f22cd..a42982e37 100644 --- a/crv/crvBezier.h +++ b/crv/crvBezier.h @@ -131,7 +131,7 @@ void getGregoryBlendedTransformationCoefficients(int blend, int type, apf::NewArray& c); /** \brief a per entity version of above */ -void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e); +void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e, bool isNew = false); /** \brief compute the matrix to transform between Bezier and Lagrange Points \details this is a support function, not actual ever needed. diff --git a/crv/crvCurveMesh.cc b/crv/crvCurveMesh.cc index 745b4b5b1..5cc39bbfb 100644 --- a/crv/crvCurveMesh.cc +++ b/crv/crvCurveMesh.cc @@ -44,7 +44,7 @@ void convertInterpolationPoints(apf::Mesh2* m, apf::MeshEntity* e, apf::destroyElement(elem); } -void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e) +void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e, bool isNew) { PCU_ALWAYS_ASSERT(m->canSnap()); int type = m->getType(e); @@ -56,9 +56,18 @@ void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e) m->setPoint(e,0,pt); return; } + // e is an edge or a face + // either way, get a length-scale by computing + // the distance b/w first two downward verts + apf::MeshEntity* down[12]; + m->getDownward(e, 0, down); + apf::Vector3 p0, p1; + m->getPoint(down[0], 0, p0); + m->getPoint(down[1], 0, p1); + double lengthScale = (p1 - p0).getLength(); apf::FieldShape * fs = m->getShape(); int non = fs->countNodesOn(type); - apf::Vector3 p, xi, pt(0,0,0); + apf::Vector3 p, xi, pt0, pt(0,0,0); for(int i = 0; i < non; ++i){ apf::ModelEntity* g = m->toModel(e); fs->getNodeXi(type,i,xi); @@ -67,7 +76,13 @@ void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e) else transferParametricOnTriSplit(m,e,xi,p); m->snapToModel(g,p,pt); - m->setPoint(e,i,pt); + if (isNew || !m->canGetClosestPoint()) { + m->setPoint(e,i,pt); + continue; + } + m->getPoint(e,i,pt0); + if (!m->isOnModel(g, pt0, lengthScale)) + m->setPoint(e,i,pt); } } diff --git a/crv/crvShapeHandler.cc b/crv/crvShapeHandler.cc index 75651f51e..9ec487ec0 100644 --- a/crv/crvShapeHandler.cc +++ b/crv/crvShapeHandler.cc @@ -149,7 +149,7 @@ class BezierTransfer : public ma::SolutionTransfer for (int i = 0; i < ne; ++i){ if(shouldSnap && midEdgeVerts[i] && isBoundaryEntity(mesh,midEdgeVerts[i])){ - snapToInterpolate(mesh,midEdgeVerts[i]); + snapToInterpolate(mesh,midEdgeVerts[i],true); } } int np = mesh->getShape()->getEntityShape(parentType)->countNodes(); @@ -173,7 +173,7 @@ class BezierTransfer : public ma::SolutionTransfer // do snapping here, inside refinement if (isBdryEnt && shouldSnap){ - snapToInterpolate(mesh,newEntities[i]); + snapToInterpolate(mesh,newEntities[i],true); } else if (useLinear && !isBdryEnt) { // boundary entities that don't get snapped should interpolate if(childType == apf::Mesh::EDGE){ @@ -538,7 +538,7 @@ class BezierHandler : public ma::ShapeHandler if (newType != apf::Mesh::EDGE) { if (snap && P > 2){ - snapToInterpolate(mesh,newEntities[i]); + snapToInterpolate(mesh,newEntities[i],true); convertInterpolationPoints(mesh,newEntities[i],n,ni,bt->coeffs[2]); } else { for (int j = 0; j < ni; ++j){ @@ -548,7 +548,7 @@ class BezierHandler : public ma::ShapeHandler } } else { if (snap){ - snapToInterpolate(mesh,newEntities[i]); + snapToInterpolate(mesh,newEntities[i],true); convertInterpolationPoints(mesh,newEntities[i],n,ni,bt->coeffs[1]); } else { setLinearEdgePoints(mesh,newEntities[i]); From 5bfcef9b61cf34a2537d0981359db500aa33c3a0 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Mon, 13 Aug 2018 15:32:05 -0400 Subject: [PATCH 18/19] fix enum for compile with gcc/4.9.2 GCC 4.9.2 does not have scoped enums, so in order to compile changes to apfConvert.cc the TagType scope needed to be removed. --- apf/apfConvert.cc | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/apf/apfConvert.cc b/apf/apfConvert.cc index ca617e06e..5b9ada911 100644 --- a/apf/apfConvert.cc +++ b/apf/apfConvert.cc @@ -8,6 +8,7 @@ #include #include #include +#include namespace apf { @@ -231,17 +232,17 @@ class Converter int* intData; long* lngData; switch (tagType) { - case apf::Mesh::TagType::DOUBLE: + case apf::Mesh::DOUBLE: dblData = new double[tagSize]; inMesh->getDoubleTag(e, in, dblData); outMesh->setDoubleTag(newFromOld[e], out, dblData); break; - case apf::Mesh::TagType::INT: + case apf::Mesh::INT: intData = new int[tagSize]; inMesh->getIntTag(e, in, intData); outMesh->setIntTag(newFromOld[e], out, intData); break; - case apf::Mesh::TagType::LONG: + case apf::Mesh::LONG: lngData = new long[tagSize]; inMesh->getLongTag(e, in, lngData); outMesh->setLongTag(newFromOld[e], out, lngData); @@ -319,13 +320,13 @@ class Converter if (!outMesh->findTag(tagName)) { apf::MeshTag* out = NULL; switch (tagType) { - case apf::Mesh::TagType::DOUBLE: + case apf::Mesh::DOUBLE: out = outMesh->createDoubleTag(tagName, tagSize); break; - case apf::Mesh::TagType::INT: + case apf::Mesh::INT: out = outMesh->createIntTag(tagName, tagSize); break; - case apf::Mesh::TagType::LONG: + case apf::Mesh::LONG: out = outMesh->createLongTag(tagName, tagSize); break; default: From 209290c919decfaacf1dc327aeedb02b3092e2df Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Mon, 13 Aug 2018 15:46:24 -0400 Subject: [PATCH 19/19] fix test case --- test/verify_convert.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/verify_convert.cc b/test/verify_convert.cc index 73d7e86e1..0ce1bb5cf 100644 --- a/test/verify_convert.cc +++ b/test/verify_convert.cc @@ -24,7 +24,7 @@ apf::Mesh2* createMesh() m->createVertex(NULL, apf::Vector3(0, 0, 0), apf::Vector3(0, 0, 0)); verts[1] = m->createVertex(NULL, apf::Vector3(1, 0, 0), apf::Vector3(1, 0, 0)); - m->createEntity(apf::Mesh::Type::EDGE, NULL, verts); + m->createEntity(apf::Mesh::EDGE, NULL, verts); return m; } class twox : public apf::Function {