diff --git a/CMakeLists.txt b/CMakeLists.txt index 70fba0900..831af6895 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -197,7 +197,6 @@ scorec_export_library(core) if(BUILD_EXES) add_subdirectory(test) - add_subdirectory(capstone_clis) endif() if(PUMI_PYTHON_INTERFACE) diff --git a/capstone_clis/CMakeLists.txt b/capstone_clis/CMakeLists.txt deleted file mode 100644 index 5d1c29f00..000000000 --- a/capstone_clis/CMakeLists.txt +++ /dev/null @@ -1,28 +0,0 @@ -function(util_exe_func exename srcname) - add_executable(${exename} ${srcname}) - target_link_libraries(${exename} core apf_cap) - bob_export_target(${exename}) - bob_end_subdir() -endfunction(util_exe_func) - -function(test_exe_func exename srcname) - if(IS_TESTING) - add_executable(${exename} ${srcname}) - else() - add_executable(${exename} EXCLUDE_FROM_ALL ${srcname}) - endif() - target_link_libraries(${exename} core) -endfunction(test_exe_func) - -# Capstone utilities -if(ENABLE_CAPSTONE) - util_exe_func(capStoneAnisoAdapt capStoneAnisoAdapt.cc) - util_exe_func(capStoneAnisoAdaptWing capStoneAnisoAdaptWing.cc) - util_exe_func(capStoneIsoAdaptB737 capStoneIsoAdaptB737.cc) - util_exe_func(capStoneGeomTest capStoneGeomTest.cc) - util_exe_func(capStoneAttachSolution capStoneAttachSolution.cc) - util_exe_func(capStoneCheckParametrization capStoneCheckParametrization.cc) - util_exe_func(capStone2VTK capStone2VTK.cc) -endif() - -bob_end_subdir() diff --git a/capstone_clis/capStoneAnalyzeArea.cc b/capstone_clis/capStoneAnalyzeArea.cc deleted file mode 100644 index fdf4a0c90..000000000 --- a/capstone_clis/capStoneAnalyzeArea.cc +++ /dev/null @@ -1,247 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - - -double totalArea(apf::Mesh2* m); - -double addSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea); - - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - lion_set_verbosity(1); - double initialTime = PCU_Time(); - - if (argc != 3) { - if(0==PCU_Comm_Self()) - std::cerr << "usage: " << argv[0] - << " \n"; - return EXIT_FAILURE; - } - - - gmi_register_mesh(); - gmi_register_null(); - - const char* creFileName = argv[1]; - int factor = atoi(argv[2]); - - // load capstone mesh - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */ - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(creFileName); - - M_GModel gmodel = cs.load_files(filenames); - - int numbreps = 0; - MG_CALL(g->get_num_breps(numbreps)); - std::cout << "number of b reps is " << numbreps << std::endl; - if(numbreps == 0) - error(HERE, ERR_INVALID_INPUT, "Model is empty"); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - for(std::size_t i = 0; i < mmodels.size(); ++i) - { - M_MModel ammodel = mmodels[i]; - std::size_t numregs = 0; - std::size_t numfaces = 0; - std::size_t numedges = 0; - std::size_t numverts = 0; - MG_API_CALL(m, set_current_model(ammodel)); - MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs)); - MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces)); - MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges)); - MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts)); - std::cout << "num regions is " << numregs << std::endl; - std::cout << "num faces is " << numfaces << std::endl; - std::cout << "num edges is " << numedges << std::endl; - std::cout << "num verts is " << numverts << std::endl; - std::cout << "-----------" << std::endl; - if(numregs > 0) - { - mmodel = ammodel; - break; - } - } - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - - gmi_cap_start(); - gmi_register_cap(); - - printf("---- CapStone Mesh Loaded. \n"); - - // create the mesh object (this one is CapStone underneath) - printf("\n---- Creating Mesh Wrapper Object. \n"); - apf::Mesh2* mesh = apf::createMesh(m,g); - printf("---- Creating Mesh Wrapper Object: Done. \n"); - - // add size distribution based on area - std::vector binCount; - std::vector binArea; - double minArea = addSizeDistribution(mesh, factor, binCount, binArea); - - printf("min area is %f\n", minArea); - - int a_lower_factor = 1; - for (std::size_t i = 0; i < binCount.size(); i++) { - double mn = i == 0 ? 0 : a_lower_factor * minArea; - double mx = i == 0 ? minArea : 4 * a_lower_factor * minArea; - printf("bin %lu's [%f,%f] count/area is %d/%g \n", i, mn, mx, binCount[i], binArea[i]); - a_lower_factor *= i == 0 ? 1 : 4; - } - - double totalBinArea = 0.0; - int totalBinCount = 0; - - for (std::size_t i = 0; i < binCount.size(); i++) { - totalBinCount += binCount[i]; - totalBinArea += binArea[i]; - } - - printf("total area is %E \n", totalArea(mesh)); - printf("total bin area is %E \n", totalBinArea); - printf("total bin count is %lu \n", totalBinCount); - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} - -double totalArea(apf::Mesh2* m) -{ - double area = 0.0; - apf::MeshEntity* e; - apf::MeshIterator* it; - it = m->begin(2); - while ( (e = m->iterate(it)) ) - area += apf::measure(m, e); - m->end(it); - return area; -} - -double addSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea) -{ - // find the min length edge - apf::MeshEntity* e; - apf::MeshIterator* it; - - double minLength = 1.0e12; - it = m->begin(1); - while ( (e = m->iterate(it)) ) { - double l = apf::measure(m, e); - if (l < minLength) - minLength = l; - } - m->end(it); - - double minArea = minLength * minLength / 2.; - - // initialize the arrays - for (int i = 1; i <= factor*factor; i*=4) { - binCount.push_back(0); - binArea.push_back(0.0); - } - - binCount.push_back(0); - binCount.push_back(0); - binArea.push_back(0.0); - binArea.push_back(0.0); - - it = m->begin(2); - while ( (e = m->iterate(it)) ) { - double a = apf::measure(m, e); - if (a < minArea){ - binCount[0] += 1; - binArea[0] += a; - continue; - } - int count = 1; - int i; - for (i = 1; i <= factor*factor; i*=4) { - if ((a >= i * minArea) && (a < 4 * i * minArea)) { - binCount[count] += 1; - binArea[count] += a; - } - count++; - } - i = i/4; - if (a >= 4 * i * minArea) { - binCount[count] += 1; - binArea[count] += a; - } - } - m->end(it); - return minArea; -} diff --git a/capstone_clis/capStoneAnalyzeEdgeLengths.cc b/capstone_clis/capStoneAnalyzeEdgeLengths.cc deleted file mode 100644 index 31dd09403..000000000 --- a/capstone_clis/capStoneAnalyzeEdgeLengths.cc +++ /dev/null @@ -1,253 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - - -double totalArea(apf::Mesh2* m); - -double addSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea); - -void edgeLengthStats(apf::Mesh2* m); - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - lion_set_verbosity(1); - double initialTime = PCU_Time(); - - if (argc != 3) { - if(0==PCU_Comm_Self()) - std::cerr << "usage: " << argv[0] - << " \n"; - return EXIT_FAILURE; - } - - - gmi_register_mesh(); - gmi_register_null(); - - const char* creFileName = argv[1]; - int factor = atoi(argv[2]); - - // load capstone mesh - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */ - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(creFileName); - - M_GModel gmodel = cs.load_files(filenames); - - int numbreps = 0; - MG_CALL(g->get_num_breps(numbreps)); - std::cout << "number of b reps is " << numbreps << std::endl; - if(numbreps == 0) - error(HERE, ERR_INVALID_INPUT, "Model is empty"); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - for(std::size_t i = 0; i < mmodels.size(); ++i) - { - M_MModel ammodel = mmodels[i]; - std::size_t numregs = 0; - std::size_t numfaces = 0; - std::size_t numedges = 0; - std::size_t numverts = 0; - MG_API_CALL(m, set_current_model(ammodel)); - MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs)); - MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces)); - MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges)); - MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts)); - std::cout << "num regions is " << numregs << std::endl; - std::cout << "num faces is " << numfaces << std::endl; - std::cout << "num edges is " << numedges << std::endl; - std::cout << "num verts is " << numverts << std::endl; - std::cout << "-----------" << std::endl; - if(numregs > 0) - { - mmodel = ammodel; - break; - } - } - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - - gmi_cap_start(); - gmi_register_cap(); - - printf("---- CapStone Mesh Loaded. \n"); - - // create the mesh object (this one is CapStone underneath) - printf("\n---- Creating Mesh Wrapper Object. \n"); - apf::Mesh2* mesh = apf::createMesh(m,g); - printf("---- Creating Mesh Wrapper Object: Done. \n"); - - // add size distribution based on area - std::vector binCount; - std::vector binArea; - edgeLengthStats(mesh); - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} - -double totalArea(apf::Mesh2* m) -{ - double area = 0.0; - apf::MeshEntity* e; - apf::MeshIterator* it; - it = m->begin(2); - while ( (e = m->iterate(it)) ) - area += apf::measure(m, e); - m->end(it); - return area; -} - -double addSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea) -{ - // find the min length edge - apf::MeshEntity* e; - apf::MeshIterator* it; - - double minLength = 1.0e12; - it = m->begin(1); - while ( (e = m->iterate(it)) ) { - double l = apf::measure(m, e); - if (l < minLength) - minLength = l; - } - m->end(it); - - double minArea = minLength * minLength / 2.; - - // initialize the arrays - for (int i = 1; i <= factor*factor; i*=4) { - binCount.push_back(0); - binArea.push_back(0.0); - } - - binCount.push_back(0); - binCount.push_back(0); - binArea.push_back(0.0); - binArea.push_back(0.0); - - it = m->begin(2); - while ( (e = m->iterate(it)) ) { - double a = apf::measure(m, e); - if (a < minArea){ - binCount[0] += 1; - binArea[0] += a; - continue; - } - int count = 1; - int i; - for (i = 1; i <= factor*factor; i*=4) { - if ((a >= i * minArea) && (a < 4 * i * minArea)) { - binCount[count] += 1; - binArea[count] += a; - } - count++; - } - i = i/4; - if (a >= 4 * i * minArea) { - binCount[count] += 1; - binArea[count] += a; - } - } - m->end(it); - return minArea; -} - -void edgeLengthStats(apf::Mesh2* m) -{ - // find the min length edge - apf::MeshEntity* e; - apf::MeshIterator* it; - - double minLength = 1.0e12; - double maxLength = 0.0; - double avgLength = 0.0; - double count = 0; - it = m->begin(1); - while ( (e = m->iterate(it)) ) { - double l = apf::measure(m, e); - if (l < minLength) - minLength = l; - if (l > maxLength) - maxLength = l; - avgLength += l; - count++; - } - m->end(it); - avgLength /= count; - - printf("min/max/avg lengths are: %E/%E/%E\n", minLength, maxLength, avgLength); - printf("num tris is %d\n", m->count(2)); -} diff --git a/capstone_clis/capStoneAnisoAdapt.cc b/capstone_clis/capStoneAnisoAdapt.cc deleted file mode 100644 index c6d71cff9..000000000 --- a/capstone_clis/capStoneAnisoAdapt.cc +++ /dev/null @@ -1,202 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - -#include "capStoneSizeFields.h" - - -void writeCre(CapstoneModule& cs, const std::string& vtkFileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, vtkFileName.c_str(), idmapping); -} - -int main(int argc, char** argv) -{ - /* INIT CALLS */ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - - if (argc < 3) { - if (PCU_Comm_Self() == 0) { - printf("USAGE: %s \n", argv[0]); - printf("Size-fields:\n"); - printf("%d, for uniform anisotropic size-field\n", 1); - printf("%d, for wing-shock size-field\n", 2); - printf("%d, for cube-shock size-field\n", 3); - } - MPI_Finalize(); - exit(EXIT_FAILURE); - } - - const char* createFileName = argv[1]; - int mode = atoi(argv[2]); - - gmi_cap_start(); - gmi_register_cap(); - - /* LOAD CAPSTONE MESH */ - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("test", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(createFileName); - - M_GModel gmodel = cs.load_files(filenames); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - PCU_ALWAYS_ASSERT(mmodels.size() == 1); - MG_API_CALL(m, set_current_model(mmodels[0])); - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - /* CONVERT THE MESH TO APF::MESH2 */ - apf::Mesh2* apfCapMesh = apf::createMesh(m, g); - - /* ADD A TEST FIELD TO THE MESH TO DEMONSTRATE SOLUTION TRANSFER */ - apf::Field* tf = apf::createFieldOn(apfCapMesh, "test_field", apf::VECTOR); - apf::MeshEntity* ent; - apf::MeshIterator* it = apfCapMesh->begin(0); - while( (ent = apfCapMesh->iterate(it)) ) { - apf::Vector3 p = ma::getPosition(apfCapMesh, ent); - double x = p[0]; - double y = p[1]; - double z = p[2]; - apf::Vector3 s(y, z, 2*x); - apf::setVector(tf, ent, 0, s); - } - apfCapMesh->end(it); - - /* WRITE THE BEFORE ADAPT MESH TO VTK USING APF VTK WRITER */ - apf::writeVtkFiles("before", apfCapMesh); - - - /* SETUP AND CALL ADAPT */ - ma::Input* in; - - // choose a size field here - ma::AnisotropicFunction* sf = 0; - switch (mode) { - case 1: - sf = new UniformAniso(apfCapMesh); - break; - case 2: - sf = new WingShock(apfCapMesh, 50); - break; - case 3: - sf = new Shock(apfCapMesh); - break; - } - - // make pumi fields that hold the "frames" and "scales" for anisotropic size fields - // here we are using user-defined size-fields. Usually, "frames" and "scales" come - // from a solution driven error estimation procedure - apf::Field* frameField = apf::createFieldOn(apfCapMesh, "adapt_frames", apf::MATRIX); - apf::Field* scaleField = apf::createFieldOn(apfCapMesh, "adapt_scales", apf::VECTOR); - - apf::MeshEntity* v; - it = apfCapMesh->begin(0); - while( (v = apfCapMesh->iterate(it)) ) { - apf::Vector3 s; - apf::Matrix3x3 f; - sf->getValue(v, f, s); - apf::setVector(scaleField, v, 0, s); - apf::setMatrix(frameField, v, 0, f); - } - apfCapMesh->end(it); - - in = ma::makeAdvanced(ma::configure(apfCapMesh, scaleField, frameField)); - in->shouldSnap = true; - in->shouldTransferParametric = true; - in->shouldFixShape = true; - in->shouldForceAdaptation = true; - if (apfCapMesh->getDimension() == 2) - in->goodQuality = 0.04; // this is mean-ratio squared - else // 3D meshes - in->goodQuality = 0.027; // this is the mean-ratio cubed - in->maximumIterations = 10; - - ma::adaptVerbose(in, false); - - /* apfCapMesh->verify(); */ - - /* WRITE THE AFTER ADAPT MESH TO VTK USING APF VTK WRITER */ - apf::writeVtkFiles("after", apfCapMesh); - - - /* PRINT ADAPTED MESH INFO */ - M_MModel mmdl; - m->get_current_model(mmdl); - std::string info; - m->print_info(mmdl, info); - - /* WRITE THE ADAPTED MESH IN NATIVE CREATE FORMAT */ - writeCre(cs, "after.cre"); - - /* CLEAN UPS */ - delete sf; - - /* EXIT CALLS */ - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} diff --git a/capstone_clis/capStoneAnisoAdaptWing.cc b/capstone_clis/capStoneAnisoAdaptWing.cc deleted file mode 100644 index 5db6666be..000000000 --- a/capstone_clis/capStoneAnisoAdaptWing.cc +++ /dev/null @@ -1,228 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" - -#include "capStoneSizeFields.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - -/* class WingShock : public ma::AnisotropicFunction */ -/* { */ -/* public: */ -/* WingShock(ma::Mesh* m, double inFactor) */ -/* { */ -/* mesh = m; */ -/* factor = inFactor; */ -/* } */ -/* virtual void getValue(ma::Entity* v, ma::Matrix& R, ma::Vector& H) */ -/* { */ -/* ma::Vector p = ma::getPosition(mesh,v); */ -/* double x = p[0]; */ -/* double x0 = 0.5; */ -/* double planeSize = 0.03125; */ -/* double spanSize = 0.5; */ -/* double delta = 0.5; */ - -/* double beta = 0.3; */ -/* double x1 = 0.2; */ -/* double f = beta + x * (1. - beta) / x1; */ -/* double multipier = (x <= x1) ? 1.0 : f; */ - -/* double s0 = planeSize / factor; */ -/* double alpha = planeSize * (1. - 1./factor) / delta; */ -/* double hx = s0 + alpha * std::abs(x - x0); */ -/* double hy = spanSize; */ -/* double hz = spanSize; */ -/* hx *= multipier; */ -/* hy *= multipier; */ - -/* ma::Vector h; */ - -/* h = ma::Vector(hx, hy, hz); */ -/* /1* // principal directions *1/ */ -/* ma::Matrix r(1.0, 0.0, 0.0, */ -/* 0.0, 1.0, 0.0, */ -/* 0.0, 0.0, 1.0); */ -/* H = h; */ -/* R = r; */ -/* } */ -/* private: */ -/* ma::Mesh* mesh; */ -/* double factor; */ -/* }; */ - -void writeCre(CapstoneModule& cs, const std::string& vtkFileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, vtkFileName.c_str(), idmapping); -} - -int main(int argc, char** argv) -{ - /* INIT CALLS */ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - - if (argc < 2) { - if (PCU_Comm_Self() == 0) { - printf("USAGE: %s \n", argv[0]); - } - MPI_Finalize(); - exit(EXIT_FAILURE); - } - - const char* createFileName = argv[1]; - - gmi_cap_start(); - gmi_register_cap(); - - /* LOAD CAPSTONE MESH */ - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("test", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(createFileName); - - M_GModel gmodel = cs.load_files(filenames); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - PCU_ALWAYS_ASSERT(mmodels.size() == 1); - MG_API_CALL(m, set_current_model(mmodels[0])); - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - /* CONVERT THE MESH TO APF::MESH2 */ - apf::Mesh2* apfCapMesh = apf::createMesh(m, g); - - /* ADD A TEST FIELD TO THE MESH TO DEMONSTRATE SOLUTION TRANSFER */ - apf::Field* tf = apf::createFieldOn(apfCapMesh, "test_field", apf::VECTOR); - apf::MeshEntity* ent; - apf::MeshIterator* it = apfCapMesh->begin(0); - while( (ent = apfCapMesh->iterate(it)) ) { - apf::Vector3 p = ma::getPosition(apfCapMesh, ent); - double x = p[0]; - double y = p[1]; - double z = p[2]; - apf::Vector3 s(y, z, 2*x); - apf::setVector(tf, ent, 0, s); - } - apfCapMesh->end(it); - - /* WRITE THE BEFORE ADAPT MESH TO VTK USING APF VTK WRITER */ - apf::writeVtkFiles("before", apfCapMesh); - - - /* SETUP AND CALL ADAPT */ - - // setting the aniso size field "scales" and "frames" - ma::IsotropicFunction* sf = 0; - sf = new Uniform(apfCapMesh, 50); - // make pumi fields that hold the "frames" and "scales" for anisotropic size fields - // here we are using user-defined size-fields. Usually, "frames" and "scales" come - // from a solution driven error estimation procedure - /* apf::Field* frameField = apf::createFieldOn(apfCapMesh, "adapt_frames", apf::MATRIX); */ - /* apf::Field* scaleField = apf::createFieldOn(apfCapMesh, "adapt_scales", apf::VECTOR); */ - apf::Field* scaleField = apf::createFieldOn(apfCapMesh, "adapt_scales", apf::SCALAR); - apf::MeshEntity* v; - it = apfCapMesh->begin(0); - while( (v = apfCapMesh->iterate(it)) ) { - /* apf::Vector3 s; */ - /* apf::Matrix3x3 f; */ - apf::setScalar(scaleField, v, 0, sf->getValue(v)); - } - apfCapMesh->end(it); - - // adapt setup - ma::Input* in; - in = ma::makeAdvanced(ma::configure(apfCapMesh, scaleField)); - /* in = ma::configure(apfCapMesh, scaleField, frameField); */ - in->shouldSnap = true; - in->shouldTransferParametric = true; - in->shouldFixShape = true; - in->shouldForceAdaptation = true; - if (apfCapMesh->getDimension() == 2) - in->goodQuality = 0.04; // this is mean-ratio squared - else // 3D meshes - in->goodQuality = 0.027; // this is the mean-ratio cubed - in->maximumIterations = 10; - - ma::adaptVerbose(in, false); - - /* WRITE THE AFTER ADAPT MESH TO VTK USING APF VTK WRITER */ - apf::writeVtkFiles("after", apfCapMesh); - - - /* PRINT ADAPTED MESH INFO */ - M_MModel mmdl; - m->get_current_model(mmdl); - std::string info; - m->print_info(mmdl, info); - - /* WRITE THE ADAPTED MESH IN NATIVE CREATE FORMAT */ - writeCre(cs, "after.cre"); - - /* CLEAN UPS */ - delete sf; - - /* EXIT CALLS */ - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} diff --git a/capstone_clis/capStoneAttachSolution.cc b/capstone_clis/capStoneAttachSolution.cc deleted file mode 100644 index 28b850a78..000000000 --- a/capstone_clis/capStoneAttachSolution.cc +++ /dev/null @@ -1,1426 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - - -typedef std::vector row; - -struct FieldOfInterest{ - FieldOfInterest() - { - lambda_max = 0.0; - h_lambdamax = 0.0; - whichLayerHasLambdaMax = NULL; - } - double lambda_max; - double h_lambdamax; - apf::Field* lambdaMaxField; - apf::Field* lambdaStrandMax; - apf::Field* sizeField; - apf::Field* whichLayerHasLambdaMax; - double lambda_cutoff() - { - return lambda_max*1e-10; - } -}; - -std::vector readTable(const char* name); - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col); - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, - std::vector > &surfToStrandMap); - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize); - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize); - -void writeCre(CapstoneModule& cs, const std::string& fileName); - -void writeMdsMesh(apf::Mesh2* m, const char* name, const char* fieldName); - -void computeSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea); - - -void adjustRefinementLevel(apf::Mesh2* m, apf::Field* finalSize, - apf::Field* currentSize, int maxLevel); - -struct SortingStruct -{ - apf::Vector3 v; - double wm; - bool operator<(const SortingStruct &other) const - { - return wm < other.wm; - } -}; - -//gradation routines from Proteus - -int gradeSizeModify(apf::Mesh* m, apf::Field* size_iso,double gradingFactor, - double size[2], apf::Adjacent edgAdjVert, - apf::Adjacent vertAdjEdg, - std::queue &markedEdges, - apf::MeshTag* isMarked, - int fieldType, - int vecPos, //which idx of sizeVec to modify - int idxFlag) - -//General function to actually modify sizes -{ - (void)vecPos; - //Determine a switching scheme depending on which vertex needs a modification - int idx1,idx2; - if(idxFlag == 0){ - idx1=0; - idx2=1; - } - else{ - idx1=1; - idx2 = 0; - } - - int marker[3] = {0,1,0}; - double marginVal = 0.01; - int needsParallel=0; - - if(fieldType == apf::SCALAR){ - //apf::Field* size_iso = m->findField("size"); - - if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal)) - { - if(m->isOwned(edgAdjVert[idx1])) - { - size[idx1] = gradingFactor*size[idx2]; - apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]); - m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - //if edge is not already marked - if(!marker[2]){ - m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]); - markedEdges.push(vertAdjEdg[i]); - } - } - } //end isOwned - else - { //Pack information to owning processor - needsParallel=1; - apf::Copies remotes; - m->getRemotes(edgAdjVert[idx1],remotes); - double newSize = gradingFactor*size[idx2]; - int owningPart=m->getOwner(edgAdjVert[idx1]); - PCU_COMM_PACK(owningPart, remotes[owningPart]); - PCU_COMM_PACK(owningPart,newSize); - } - } - - }//end if apf::SCALAR - return needsParallel; -} - -void markEdgesInitial(apf::Mesh* m, apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used to initially determine which edges need to be considered for gradation -{ - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - double size[2]; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::MeshEntity* edge; - apf::MeshIterator* it = m->begin(1); - while((edge=m->iterate(it))){ - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i]=apf::getScalar(size_iso,edgAdjVert[i],0); - } - if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ - //add edge to a queue - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - else{ - m->setIntTag(edge,isMarked,&marker[0]); - } - } - m->end(it); -} - -int serialGradation(apf::Mesh* m,apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used loop over the mesh edge queue for gradation and modify the sizes -{ - double size[2]; - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - apf::MeshEntity* edge; - int needsParallel=0; - - //perform serial gradation while packing necessary info for parallel - while(!markedEdges.empty()){ - edge = markedEdges.front(); - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i] = apf::getScalar(size_iso,edgAdjVert[i],0); - } - - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 0); - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 1); - - m->setIntTag(edge,isMarked,&marker[0]); - markedEdges.pop(); - } - return needsParallel; -} - -int gradeMesh(apf::Mesh* m,apf::Field* size_iso) -//Function to grade isotropic mesh through comparison of edge vertex size ratios -//This implementation accounts for parallel meshes as well -//First do serial gradation. -//If a shared entity has its size modified, then send new size to owning copy. -//After full loop over entities, have owning copy take minimum of all sizes received -//Flag adjacent entities to owning copy. -//Communicate to remote copies that a size was modified, and flag adjacent edges to remote copies for further gradation -{ - //unique to this code - //apf::Field* size_iso = m->findField("size"); - double gradingFactor = 1.2; - // - - apf::MeshEntity* edge; - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - std::queue markedEdges; - apf::MeshTag* isMarked = m->createIntTag("isMarked",1); - - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - apf::MeshIterator* it; - markEdgesInitial(m,size_iso,markedEdges,gradingFactor); - - int needsParallel=1; - while(needsParallel) - { - PCU_Comm_Begin(); - needsParallel = serialGradation(m,size_iso,markedEdges,gradingFactor); - - PCU_Add_Ints(&needsParallel,1); - PCU_Comm_Send(); - - apf::MeshEntity* ent; - double receivedSize; - double currentSize; - double newSize; - - //Need a container to get all entitites that need to be updated on remotes - std::queue updateRemoteVertices; - - apf::Copies remotes; - //owning copies are receiving - while(PCU_Comm_Receive()) - { - PCU_COMM_UNPACK(ent); - PCU_COMM_UNPACK(receivedSize); - - if(!m->isOwned(ent)){ - std::cout<<"THERE WAS AN ERROR"<getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - updateRemoteVertices.push(ent); - } - - PCU_Comm_Begin(); - - while(!updateRemoteVertices.empty()) - { - ent = updateRemoteVertices.front(); - //get remote copies and send updated mesh sizes - m->getRemotes(ent,remotes); - currentSize = apf::getScalar(size_iso,ent,0); - for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter) - { - PCU_COMM_PACK(iter->first, iter->second); - } - updateRemoteVertices.pop(); - } - - PCU_Comm_Send(); - //while remote copies are receiving - while(PCU_Comm_Receive()) - { - //unpack - PCU_COMM_UNPACK(ent); - //PCU_COMM_UNPACK(receivedSize); - assert(!m->isOwned(ent)); - - if(m->isOwned(ent)){ - std::cout<<"Problem occurred\n"; - std::exit(1); - } - - //add adjacent edges into Q - m->getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - } - apf::synchronize(size_iso); - - } //end outer while - - //Cleanup of edge marker field - it = m->begin(1); - while((edge=m->iterate(it))){ - m->removeTag(edge,isMarked); - } - m->end(it); - m->destroyTag(isMarked); - - //apf::synchronize(size_iso); - return needsParallel; -} - -//Eigenvalue routines for volume and surface meshes -// - //function to get volume-based eigenvalues -//apf::Field* getLambdaMax(mesh,hessianField,) -void getLambdaMax(apf::Mesh* mesh,apf::Field* hessianField,apf::Field* lambdaMaxField, apf::Field* lambdaMaxEVecField = NULL) -{ - apf::MeshEntity* vert; - apf::MeshIterator* it; - it = mesh->begin(0); - //apf::Field* lambdaMaxField = apf::createLagrangeField(mesh,"lambdaMax",apf::SCALAR,1); - while( (vert = mesh->iterate(it)) ) - { - apf::Matrix3x3 metric; - apf::getMatrix(hessianField, vert, 0, metric); - apf::Vector3 eigenVectors[3]; - double eigenValues[3]; - metric = (metric+apf::transpose(metric))/2.0; - apf::eigen(metric, eigenVectors, eigenValues); - - // Sort the eigenvalues and corresponding vectors - // Larger eigenvalues means a need for a finer mesh - SortingStruct ssa[3]; - for (int i = 0; i < 3; ++i) - { - ssa[i].v = eigenVectors[i]; - ssa[i].wm = std::fabs(eigenValues[i]); - } - std::sort(ssa, ssa + 3); - - assert(ssa[2].wm >= ssa[1].wm); - assert(ssa[1].wm >= ssa[0].wm); - apf::setScalar(lambdaMaxField,vert,0,ssa[2].wm); - if (lambdaMaxEVecField) - { - apf::setVector(lambdaMaxEVecField,vert,0,ssa[2].v); - } - } - mesh->end(it); -} - -void getVolMaxPair(apf::Mesh* mesh,std::vector > surfToStrandMap, apf::Field* lambdaMaxField,apf::Field* lambdaStrandMax, apf::Field* currentSize,double &lambda_max,double &h_lambdamax, apf::Field* whichLayerHasLambdaMax = NULL) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert, *vertVol; - int counter = 0; - while( (vert = mesh->iterate(it)) ) - { - //the goal of this loop is to identify lambda max global, find corresponding surface h-max global, and then set the surface lambda field - int strandSize = surfToStrandMap[counter].size(); - for(int i=0; iend(it); -} - -void getSurfMaxPair(apf::Mesh* mesh,apf::Field* lambdaMaxField,apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert; - while( (vert = mesh->iterate(it)) ) - { - if(lambda_max < apf::getScalar(lambdaMaxField,vert,0)) - { - lambda_max = apf::getScalar(lambdaMaxField,vert,0); - h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - mesh->end(it); -} - -void setSizeField(apf::Mesh* mesh, apf::Field* lambdaMaxField,apf::Field* sizeField,double lambda_max,double lambda_cutoff, double h_lambdamax,double h_global,double factor) -{ - apf::MeshIterator* it=mesh->begin(0); - apf::MeshEntity* vert; - int counter3 = 0; - double h_special = -1; - while( (vert = mesh->iterate(it)) ) - { - double h_v = h_special; //set default size as user specified input - double lambda_vert = apf::getScalar(lambdaMaxField,vert,0); - if(lambda_vert > lambda_cutoff) - { - h_v = sqrt(lambda_max/lambda_vert/factor)*h_lambdamax; - } - else - counter3++; - //h_global is user-specified max size - if(h_v > h_global) //maximum value for size field - h_v = h_global; - apf::setScalar(sizeField,vert,0,h_v); - } - mesh->end(it); - -} - -void setSizeFieldAlt(apf::Mesh* mesh, apf::Field* lambdaMaxField,apf::Field* sizeField,apf::Field* currentSize,double lambda_max,double lambda_cutoff,double h_global,double factor) -{ - apf::MeshIterator* it=mesh->begin(0); - apf::MeshEntity* vert; - int counter3 = 0; - double h_special = -1; - while( (vert = mesh->iterate(it)) ) - { - double h_v = h_special; //set default size as user specified input - double lambda_vert = apf::getScalar(lambdaMaxField,vert,0); - double h_v_curr = apf::getScalar(currentSize,vert,0); - if(lambda_vert > lambda_cutoff) - { - h_v = sqrt(lambda_max/lambda_vert/factor)*h_v_curr; - } - else - { - h_v = h_global; - } - //h_global is user-specified max size - if(h_v > h_global) //maximum value for size field - h_v = h_global; - apf::setScalar(sizeField,vert,0,h_v); - } - mesh->end(it); - -} - -void setSizeFieldAlt2(apf::Mesh* mesh, apf::Field* lambdaMaxField,apf::Field* sizeField,apf::Field* currentSize,double lambda_cutoff,double h_global,double factor) -{ - apf::MeshIterator* it=mesh->begin(0); - apf::MeshEntity* vert; - int counter3 = 0; - double h_special = -1; - while( (vert = mesh->iterate(it)) ) - { - double h_v = h_special; //set default size as user specified input - double lambda_vert = apf::getScalar(lambdaMaxField,vert,0); - double h_v_curr = apf::getScalar(currentSize,vert,0); - if(lambda_vert > lambda_cutoff) - { - h_v = h_v_curr / factor; - } - else - { - h_v = h_global; - } - //h_global is user-specified max size - if(h_v > h_global) //maximum value for size field - h_v = h_global; - apf::setScalar(sizeField,vert,0,h_v); - } - mesh->end(it); - -} - -static std::vector decodeBitFields(const char* bitFields) -{ - std::vector res; - res.resize(strlen(bitFields)); - printf("bitfields length is %zd\n", strlen(bitFields)); - PCU_ALWAYS_ASSERT_VERBOSE(strlen(bitFields) == 4, "bitfields length needs to be 4"); - for (unsigned i = 0; i < strlen(bitFields); i++) { - if (bitFields[i] == '0') - res[i] = false; - else if (bitFields[i] == '1') - res[i] = true; - else - PCU_ALWAYS_ASSERT_VERBOSE(0, "Bit fields need to be 0 or 1"); - } - PCU_ALWAYS_ASSERT_VERBOSE((bitFields[2] == '0') || (bitFields[3] == '0'), - "Both bit field 2 and bit field 3 cannot be 1"); - return res; -} - -void isotropicIntersect(apf::Mesh* m, std::queue sizeFieldList, - const std::vector& userInput, apf::Field* finalSizeField, - apf::Field* finalChoiceField) -{ - // std::vector userInput = decodeBitFields(bitFields); - PCU_ALWAYS_ASSERT(userInput.size() == sizeFieldList.size()); - apf::MeshEntity *vert; - - apf::Field *field = sizeFieldList.front(); - apf::MeshIterator *it = m->begin(0); - double largeNum = 1.e16; - while( (vert = m->iterate(it)) ) - { - apf::setScalar(finalSizeField,vert,0,largeNum); - apf::setScalar(finalChoiceField,vert,0,-1); - } - m->end(it); - int choiceIdx = 0; //assumes the initial field was set to choice 0 - while(!sizeFieldList.empty()) - { - if (userInput[choiceIdx]) { - field = sizeFieldList.front(); - apf::MeshIterator *it = m->begin(0); - while( (vert = m->iterate(it)) ) - { - double value1 = apf::getScalar(finalSizeField,vert,0); - double value2 = apf::getScalar(field,vert,0); - double minValue = std::min(value1,value2); - apf::setScalar(finalSizeField,vert,0,minValue); - if(value1 > value2) - { - apf::setScalar(finalChoiceField,vert,0,choiceIdx); - } - } - m->end(it); - } - sizeFieldList.pop(); - choiceIdx++; - } -} - - -// - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - lion_set_verbosity(1); - double initialTime = PCU_Time(); - - if (argc != 9) { - if(0==PCU_Comm_Self()) { - std::cerr << "usage: " << argv[0] - << " \n"; - std::cerr << "*target field(s) is a bit string to select which field(s) are used for error estimation\n"; - std::cerr << "1st bit --- pressure\n"; - std::cerr << "2nd bit --- skin friction\n"; - } - return EXIT_FAILURE; - } - - - gmi_register_mesh(); - gmi_register_null(); - - const char* creFileName = argv[1]; - const char* dataFileName = argv[2]; - const char* bitFields = argv[3]; - const int strandSize = atoi(argv[4]); - double h_global = atof(argv[5]); - double factor = atof(argv[6]); - const int maxLevel = atoi(argv[7]); - const int nIgnoredLayers = atoi(argv[8]); - - // Some input validation - std::vector userInputFields = decodeBitFields(bitFields); - - // load capstone mesh - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */ - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(creFileName); - - M_GModel gmodel = cs.load_files(filenames); - - int numbreps = 0; - MG_CALL(g->get_num_breps(numbreps)); - std::cout << "number of b reps is " << numbreps << std::endl; - if(numbreps == 0) - error(HERE, ERR_INVALID_INPUT, "Model is empty"); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - for(std::size_t i = 0; i < mmodels.size(); ++i) - { - M_MModel ammodel = mmodels[i]; - std::size_t numregs = 0; - std::size_t numfaces = 0; - std::size_t numedges = 0; - std::size_t numverts = 0; - MG_API_CALL(m, set_current_model(ammodel)); - MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs)); - MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces)); - MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges)); - MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts)); - std::cout << "num regions is " << numregs << std::endl; - std::cout << "num faces is " << numfaces << std::endl; - std::cout << "num edges is " << numedges << std::endl; - std::cout << "num verts is " << numverts << std::endl; - std::cout << "-----------" << std::endl; - if(numregs > 0) - { - mmodel = ammodel; - break; - } - } - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - - gmi_cap_start(); - gmi_register_cap(); - - printf("---- CapStone Mesh Loaded. \n"); - - // read the data into a vector array for now - printf("\n---- Reading Strand Data. \n"); - std::vector table = readTable(dataFileName); - printf("---- Reading Strand Data: Done. \n"); - - // create the mesh object (this one is CapStone underneath) - printf("\n---- Creating Mesh Wrapper Object. \n"); - apf::Mesh2* mesh = apf::createMesh(m,g); - printf("---- Creating Mesh Wrapper Object: Done. \n"); - - // make the volume mesh (this one is MDS underneath) - printf("\n---- Creating Volume Mesh. \n"); - // create an empty array with the correct sizes to hold the surface 2 strand map - std::vector > surfToStrandMap(mesh->count(0), - std::vector(strandSize, NULL)); - apf::Mesh2* volMesh = createVolumeMesh(mesh, table, strandSize, surfToStrandMap); - printf("---- Creating Volume Mesh: Done. \n"); - - // get all fields and add them to the mesh - printf("\n---- Adding Fields to Volume Mesh. \n"); - apf::Field* rhoField = addScalarField(volMesh, table, "rho", 3, strandSize); - apf::Field* rho_uvwField = addVector3Field(volMesh, table, "rho_uvw", 4, 5, 6, strandSize); - apf::Field* eField = addScalarField(volMesh, table, "e", 7, strandSize); - apf::Field* nuField = addScalarField(volMesh, table, "nu", 8, strandSize); - printf("---- Adding Fields to Volume Mesh: Done. \n"); - - printf("\n---- Printing Volume/Strand Mesh Stats. \n"); - printf("number of mesh verts: %zu\n", volMesh->count(0)); - printf("number of mesh regions(hexes): %zu\n", volMesh->count(3)); - printf("---- Printing Volume/Strand Mesh Stats: Done. \n"); - - double constructionTime = PCU_Time(); - std::cout<<"TIMER: Finished converting capstone mesh to volume mesh "<begin(0); - while( (vert = volMesh->iterate(it)) ) - { - double rhoVal = apf::getScalar(rhoField,vert,0); - apf::Vector3 rhoVelocity; - apf::getVector(rho_uvwField,vert,0,rhoVelocity); - double speed = rhoVelocity.getLength()/rhoVal; - apf::setScalar(speedField,vert,0,speed); - } - volMesh->end(it); - // Also get its Hessian - apf::Field* gradSpeedField = apf::recoverGradientByVolume(speedField); - apf::Field* hessianSpeedField = apf::recoverGradientByVolume(gradSpeedField); - //End getSpeed - - apf::Field* lmaxHSpeedH2Field = apf::createLagrangeField(volMesh,"lmax_hess_speed_x_h2",apf::SCALAR,1); - apf::Field* lhSpeedDotNormal = apf::createLagrangeField(volMesh,"lhSpeedDotNormal",apf::SCALAR,1); - - //get eigenvalues in the volume mesh - getLambdaMax(volMesh,hessianSpeedField,speedLambdaMaxField,speedLambdaMaxEVecField); - // NOTE: Following loop assumes that each vector in surfToStrandMap contains - // vertices starting from layer_num 1 to 51 along a strand - APF_ITERATE(std::vector >,surfToStrandMap,it) - { - double lhSpeed; - double h; - apf::Vector3 thisPoint, lastPoint, nextPoint; - apf::Vector3 diff0, diff1; - apf::Vector3 evec, strand_dir; - volMesh->getPoint((*it)[0], 0, thisPoint); - volMesh->getPoint((*it)[1], 0, nextPoint); - diff1 = nextPoint - thisPoint; - lhSpeed = apf::getScalar(speedLambdaMaxField, (*it)[0], 0); - strand_dir = diff1; - h = diff1.getLength(); - apf::setScalar(lmaxHSpeedH2Field, (*it)[0], 0, lhSpeed * h * h); - - apf::getVector(speedLambdaMaxEVecField, (*it)[0], 0, evec); - evec = evec.normalize(); - strand_dir = strand_dir.normalize(); - apf::setScalar(lhSpeedDotNormal, (*it)[0], 0, std::fabs(strand_dir * evec)); - - for (int i = 1; i < (strandSize - 1); i++) - { - lastPoint = thisPoint; - thisPoint = nextPoint; - volMesh->getPoint((*it)[i + 1], 0, nextPoint); - diff0 = thisPoint - lastPoint; - diff1 = nextPoint - thisPoint; - strand_dir = nextPoint - lastPoint; - lhSpeed = apf::getScalar(speedLambdaMaxField, (*it)[i], 0); - h = (diff0.getLength() + diff1.getLength()) / 2.0; - apf::setScalar(lmaxHSpeedH2Field, (*it)[i], 0, lhSpeed * h * h); - - apf::getVector(speedLambdaMaxEVecField, (*it)[i], 0, evec); - evec = evec.normalize(); - strand_dir = strand_dir.normalize(); - apf::setScalar(lhSpeedDotNormal, (*it)[i], 0, std::fabs(strand_dir * evec)); - } - - lastPoint = thisPoint; - thisPoint = nextPoint; - diff0 = thisPoint - lastPoint; - strand_dir = thisPoint - lastPoint; - lhSpeed = apf::getScalar(speedLambdaMaxField, (*it)[strandSize - 1], 0); - h = diff0.getLength(); - apf::setScalar(lmaxHSpeedH2Field, (*it)[strandSize - 1], 0, lhSpeed * h * h); - - apf::getVector(speedLambdaMaxEVecField, (*it)[strandSize - 1], 0, evec); - evec = evec.normalize(); - strand_dir = strand_dir.normalize(); - apf::setScalar(lhSpeedDotNormal, (*it)[strandSize - 1], 0, std::fabs(strand_dir * evec)); - } - - FieldOfInterest speedBased; - speedBased.lambdaMaxField = lmaxHSpeedH2Field; - speedBased.lambdaStrandMax = apf::createLagrangeField(mesh,"surf_lhSpeed_x_h2_strandMax",apf::SCALAR,1); - speedBased.whichLayerHasLambdaMax = apf::createLagrangeField(mesh,"surf_idx_lhSpeed_x_h2_strandMax",apf::SCALAR,1); - speedBased.sizeField = apf::createLagrangeField(mesh,"surface_size_lhSpeed_x_h2",apf::SCALAR,1); - - getVolMaxPair(mesh,surfToStrandMap,speedBased.lambdaMaxField,speedBased.lambdaStrandMax,currentSize,speedBased.lambda_max,speedBased.h_lambdamax,speedBased.whichLayerHasLambdaMax); - - // printf("before->speedBased.lambda_max: %15.10e\n", speedBased.lambda_max); - // printf("before->speedBased.h_lambdamax: %15.10e\n", speedBased.h_lambdamax); - double cutoff = speedBased.lambda_max; - speedBased.h_lambdamax = 0.0; - speedBased.lambda_max = 0.0; - it = mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - double idx_with_lmax = apf::getScalar(speedBased.whichLayerHasLambdaMax, vert, 0); - if (idx_with_lmax < nIgnoredLayers) - { - apf::setScalar(speedBased.lambdaStrandMax, vert, 0, 0.0); - } - else if (speedBased.lambda_max < apf::getScalar(speedBased.lambdaStrandMax, vert, 0)) - { - speedBased.lambda_max = apf::getScalar(speedBased.lambdaStrandMax, vert, 0); - // printf("->speedBased.lambda_max: %15.10e\n", speedBased.lambda_max); - // printf("->speedBased.lambda_max should be: %15.10e\n", apf::getScalar(speedBased.lambdaStrandMax, vert, 0)); - speedBased.h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - mesh->end(it); - // printf("after->speedBased.lambda_max: %15.10e\n", speedBased.lambda_max); - // printf("after->speedBased.h_lambdamax: %15.10e\n", speedBased.h_lambdamax); - - // Using a separate size field computation to capture free shear layer - if (bitFields[2] == '1') { - setSizeFieldAlt(mesh,speedBased.lambdaStrandMax,speedBased.sizeField,currentSize,speedBased.lambda_max,speedBased.lambda_cutoff(),h_global,factor); - } else if (bitFields[3] == '1') { - // TODO: Change from hard-coding - double factor_fsl = 4.0; - setSizeFieldAlt2(mesh,speedBased.lambdaStrandMax,speedBased.sizeField,currentSize,speedBased.lambda_cutoff(),h_global,factor_fsl); - } - - apf::Field* surfaceSpeedField = apf::createLagrangeField(mesh,"surface_speed",apf::SCALAR,1); - it = mesh->begin(0); - apf::MeshIterator* itVol = volMesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::MeshEntity* vertVol = volMesh->iterate(itVol); - double rhoVal = apf::getScalar(rhoField,vertVol,0); - double eVal = apf::getScalar(eField,vertVol,0); - double speedVal = apf::getScalar(speedField,vertVol,0); - apf::setScalar(surfaceRhoField,vert,0,rhoVal); - apf::setScalar(surfaceEField,vert,0,eVal); - apf::setScalar(surfaceSpeedField,vert,0,speedVal); - } - mesh->end(it); - volMesh->end(itVol); - std::cout<<"Finished surface fields from volume\n"; - - double getSurfaceTime = PCU_Time(); - std::cout<<"TIMER: Finished computing speed and transferring fields to surface "<begin(0); - while( (vert = volMesh->iterate(itVol)) ) - { - double rho_e = apf::getScalar(eField,vert,0); - double rho = apf::getScalar(rhoField,vert,0); - double speedVal = apf::getScalar(speedField,vert,0); - double pressure_gamma = rho_e-0.5*rho*speedVal*speedVal; - apf::setScalar(pressureField,vert,0,pressure_gamma); - } - volMesh->end(itVol); - apf::Field* gradPressureField = apf::recoverGradientByVolume(pressureField); - apf::Field* hessianPressureField = apf::recoverGradientByVolume(gradPressureField); - - //get eigenvalues in the volume mesh - getLambdaMax(volMesh,hessianPressureField,lambdaMaxField); - - - FieldOfInterest eBased; - eBased.lambdaMaxField = lambdaMaxField; - eBased.lambdaStrandMax = apf::createLagrangeField(mesh,"surf_lambda_strandMax",apf::SCALAR,1); - eBased.sizeField = apf::createLagrangeField(mesh,"surface_size_e",apf::SCALAR,1); - - getVolMaxPair(mesh,surfToStrandMap,eBased.lambdaMaxField,eBased.lambdaStrandMax,currentSize,eBased.lambda_max,eBased.h_lambdamax); - - //set size field - - setSizeField(mesh,eBased.lambdaStrandMax,eBased.sizeField,eBased.lambda_max,eBased.lambda_cutoff(),eBased.h_lambdamax,h_global,factor); - - double getPressureTime = PCU_Time(); - std::cout<<"TIMER: Finished pressure size field "<begin(0); - int vIDcounter = 0; //better way to do this would be with a numbering system - //apf::Numbering* nVID = volMesh->findNumbering; //this would be for a volume mesh - - while( (vert = mesh->iterate(it)) ) - { - apf::Vector3 pointVol; - apf::Vector3 pointSurf; - apf::MeshEntity* volVert = surfToStrandMap[vIDcounter][1]; - double speed = apf::getScalar(speedField,volVert,0); - mesh->getPoint(vert,0,pointSurf); - volMesh->getPoint(volVert,0,pointVol); - apf::Vector3 differenceVec = pointVol-pointSurf; - double shearStress = speed/differenceVec.getLength(); - apf::setScalar(surfaceShearField,vert,0,shearStress); - vIDcounter++; - } - mesh->end(it); - - FieldOfInterest shearBased; - shearBased.lambdaMaxField = apf::createLagrangeField(mesh,"surflambdaMax",apf::SCALAR,1); - shearBased.sizeField = apf::createLagrangeField(mesh,"surface_size",apf::SCALAR,1); - - //get surface shear gradient - std::cout<<"Reached gradshearfield\n"; - apf::Field* gradShearField = apf::recoverGradientByVolume(surfaceShearField); - std::cout<<"Got gradshearfield\n"; - - it=mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::Vector3 gradShear; - apf::getVector(gradShearField,vert,0,gradShear); - apf::setScalar(shearBased.lambdaMaxField,vert,0,gradShear.getLength()); - } - mesh->end(it); - - std::cout<<"Got surf lambda max\n"; - getSurfMaxPair(mesh,shearBased.lambdaMaxField,currentSize,shearBased.lambda_max,shearBased.h_lambdamax); - - std::cout<<"Got surf lambda max pair\n"; - //set size field - setSizeField(mesh,shearBased.lambdaMaxField,shearBased.sizeField,shearBased.lambda_max,shearBased.lambda_cutoff(),shearBased.h_lambdamax,h_global,factor); - - std::cout<<"set size field\n"; - - double getShearTime = PCU_Time(); - std::cout<<"TIMER: Finished skin friction size field "< sizeFieldList; - sizeFieldList.push(eBased.sizeField); - sizeFieldList.push(shearBased.sizeField); - sizeFieldList.push(speedBased.sizeField); - sizeFieldList.push(speedBased.sizeField); - apf::Field* finalChoiceField = apf::createLagrangeField(mesh,"finalChoice",apf::SCALAR,1); - it = mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::setScalar(finalChoiceField,vert,0,0); - } - - isotropicIntersect(mesh,sizeFieldList,userInputFields,finalSizeField,finalChoiceField); - - double getIntersectionTime = PCU_Time(); - std::cout<<"TIMER: get mesh intersection "<shouldRunPreZoltan = true; - in->shouldRunMidZoltan = true; - in->shouldRunPostZoltan = true; - in->maximumImbalance = 1.05; - in->maximumIterations = log2(factor) + 2; - in->shouldSnap = true; - in->shouldTransferParametric = true; - /* in->shouldTransferToClosestPoint = true; */ - in->shouldForceAdaptation = true; - /* in->debugFolder = "debug"; */ - ma::adaptVerbose(in, false); - - double adaptTime = PCU_Time(); - std::cout<<"TIMER: adaptMesh "< binCount; - std::vector binArea; - computeSizeDistribution(mesh, factor, binCount, binArea); - - for (std::size_t i = 0; i < binCount.size(); i++) { - printf("bin %lu's count/area is %d/%g \n", i, binCount[i], binArea[i]); - } - - - // write the adapted mesh to vtk - apf::writeVtkFiles("adaptedMesh", mesh); - // write the adapted mesh to smb - // Note: This is for statistic measurements using - // measureIsoStats.cc - // The last argument is the name of the field to be used - // in measureIsoStats.cc - writeMdsMesh(mesh, "adaptedMesh.smb", "adapt_size"); - // write the adapted mesh to cre - writeCre(cs, "adaptedMesh.cre"); - - // clean up and exit calls - //destroy all fields... - for(int i =0;icountFields();i++) - { - apf::destroyField(mesh->getField(i)); - } - for(int i =0;icountFields();i++) - { - apf::destroyField(volMesh->getField(i)); - } - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} - -std::vector readTable(const char* name) -{ - std::vector table; - std::ifstream file(name); - while(!file.eof()){ - std::string str; - std::getline(file, str); - std::stringstream ss(str); - double value; - row r; - while (ss >> value) - r.push_back(value); - if (r.size() > 0) - table.push_back(r); - } - return table; -} - -static std::vector extractNthStrandData(const std::vector &table, - const int strandSize, const int n, const int col) -{ - std::vector outData; - int i = 0; - while (i*strandSize+(strandSize-n) < int(table.size())) { - outData.push_back(table[i*strandSize+(strandSize - n)][col]); - i++; - } - return outData; -} - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col) -{ - // surface data corresponds to strand 1; - return extractNthStrandData(table, strandSize, 1, col); -} - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::SCALAR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value = t[vid*strandSize+(strandSize - layer)][col]; - apf::setScalar(f, e, 0, value); - count++; - } - m->end(it); - return f; -} - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::VECTOR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value0 = t[vid*strandSize+(strandSize - layer)][col0]; - double value1 = t[vid*strandSize+(strandSize - layer)][col1]; - double value2 = t[vid*strandSize+(strandSize - layer)][col2]; - apf::setVector(f, e, 0, apf::Vector3(value0, value1, value2)); - count++; - } - m->end(it); - return f; -} - - - -static std::vector -readLayerCoordinates(int layer, const std::vector &t, int s) -{ - std::vector xs = extractNthStrandData(t, s, layer, 0); - std::vector ys = extractNthStrandData(t, s, layer, 1); - std::vector zs = extractNthStrandData(t, s, layer, 2); - - std::vector res; - - for (std::size_t i = 0; i < xs.size(); i++) - res.push_back(apf::Vector3(xs[i], ys[i], zs[i])); - - return res; -} - - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, std::vector > &surfToStrandMap) -{ - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "surf_verts_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - // get the original faces - std::vector faces; - it = m->begin(2); - while ( (e = m->iterate(it)) ) - faces.push_back(e); - m->end(it); - - // create and empty mds mesh - apf::Mesh2* vMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 3, false); - apf::Numbering* nLayer = apf::createNumbering(vMesh, "layer_num", apf::getLagrange(1) , 1); - apf::Numbering* nVID = apf::createNumbering(vMesh, "vid_num", apf::getLagrange(1) , 1); - - // process each layer - // layer = 1 corresponds to the surface nodes - // treat it out side of the loop - std::vector coords = readLayerCoordinates(1, t, s); - std::vector lastVs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, 1); - apf::number(nVID, v, 0, 0, i); - lastVs.push_back(v); - } - - PCU_ALWAYS_ASSERT(coords.size() == lastVs.size()); - - // now process the remaining layers - for (int layer = 2; layer < s+1; layer++) { - // read the coords first - std::vector coords = readLayerCoordinates(layer, t, s); - // create the new verts - std::vector vs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, layer); - apf::number(nVID, v, 0, 0, i); - vs.push_back(v); - } - - // create the prisms - for (std::size_t i = 0; i < faces.size(); i++) { - apf::MeshEntity* f = faces[i]; - apf::MeshEntity* dv[3]; - m->getDownward(f, 0, dv); - int id0 = apf::getNumber(n, dv[0], 0, 0); - int id1 = apf::getNumber(n, dv[1], 0, 0); - int id2 = apf::getNumber(n, dv[2], 0, 0); - - apf::MeshEntity* prismVs[6]; - - prismVs[0] = lastVs[id0]; - prismVs[1] = lastVs[id1]; - prismVs[2] = lastVs[id2]; - - prismVs[3] = vs[id0]; - prismVs[4] = vs[id1]; - prismVs[5] = vs[id2]; - - apf::buildElement(vMesh, 0, apf::Mesh::PRISM, prismVs); - } - // set the lastVs to Vs - for (std::size_t i = 0; i < vs.size(); i++) - lastVs[i] = vs[i]; - } - vMesh->acceptChanges(); - apf::deriveMdsModel(vMesh); - apf::verify(vMesh); - - // create the surface to strand map - it = vMesh->begin(0); - while ( (e = vMesh->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0) - 1; - int vid = apf::getNumber(nVID, e, 0, 0); - PCU_ALWAYS_ASSERT((layer >= 0) && (layer < s)); - PCU_ALWAYS_ASSERT((vid >= 0) && (vid < (int)m->count(0))); - surfToStrandMap[vid][layer] = e; - } - vMesh->end(it); - - // return the volume mesh - return vMesh; -} - -void writeCre(CapstoneModule& cs, const std::string& fileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, fileName.c_str(), idmapping); -} - -void writeMdsMesh(apf::Mesh2* m, const char* name, const char* fieldName) -{ - // this would be the smb mesh with a field associated with it - apf::Mesh2* outMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 2, false); - - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "v_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - int numVerts = m->count(0); - apf::MeshEntity* vMap[numVerts]; - - apf::MeshEntity* v; - it = m->begin(0); - while( (v = m->iterate(it)) ) { - apf::Vector3 p; - m->getPoint(v, 0, p); - apf::MeshEntity* nv = outMesh->createVertex(0, p, apf::Vector3(0, 0, 0)); - int id = apf::getNumber(n, v, 0, 0); - vMap[id] = nv; - } - m->end(it); - - it = m->begin(2); - while( (e = m->iterate(it)) ) { - apf::MeshEntity* vs[3]; - m->getDownward(e, 0, vs); - apf::MeshEntity* triVs[3]; - triVs[0] = vMap[apf::getNumber(n, vs[0], 0, 0)]; - triVs[1] = vMap[apf::getNumber(n, vs[1], 0, 0)]; - triVs[2] = vMap[apf::getNumber(n, vs[2], 0, 0)]; - apf::buildElement(outMesh, 0, apf::Mesh::TRIANGLE, triVs); - } - m->end(it); - - outMesh->acceptChanges(); - apf::deriveMdsModel(outMesh); - outMesh->verify(); - apf::verify(outMesh); - - apf::destroyNumbering(n); - - // get the field on input mesh - apf::Field* inputField = m->findField(fieldName); - PCU_ALWAYS_ASSERT(inputField); - - // create the corresponding field for the out mesh - apf::Field* outputField = apf::createFieldOn(outMesh, "sizes", apf::SCALAR); - - it = m->begin(0); - apf::MeshIterator* outIt = outMesh->begin(0); - while( (v = m->iterate(it)) ) { - double value = apf::getScalar(inputField, v, 0); - apf::MeshEntity* outV = outMesh->iterate(outIt); - apf::setScalar(outputField, outV, 0, value); - } - m->end(it); - outMesh->end(outIt); - outMesh->writeNative(name); - - // clean up outMesh - outMesh->destroyNative(); - apf::destroyMesh(outMesh); -} - - -void computeSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea) -{ - // find the min length edge - apf::MeshEntity* e; - apf::MeshIterator* it; - - double minLength = 1.0e12; - it = m->begin(1); - while ( (e = m->iterate(it)) ) { - double l = apf::measure(m, e); - if (l < minLength) - minLength = l; - } - m->end(it); - - double minArea = minLength * minLength / 2.; - - // initialize the arrays - for (int i = 1; i <= factor; i*=2) { - binCount.push_back(0); - binArea.push_back(0.0); - } - - - it = m->begin(2); - while ( (e = m->iterate(it)) ) { - double a = apf::measure(m, e); - int count = 0; - for (int i = 1; i <= factor; i*=2) { - if (a >= i * minArea && a < 4 * i * minArea) { - binCount[count] += 1; - binArea[count] += a; - break; - } - count++; - } - } - m->end(it); -} - -void adjustRefinementLevel(apf::Mesh2* m, apf::Field* finalSize, - apf::Field* currentSize, int maxLevel) -{ - apf::MeshEntity* v; - apf::MeshIterator* it = m->begin(0); - - while( (v = m->iterate(it)) ) { - double currnetS = apf::getScalar(currentSize, v, 0); - double finalS = apf::getScalar(finalSize, v, 0); - if (currnetS < finalS) continue; - if (log2(currnetS/finalS) > maxLevel) - finalS = currnetS / pow(2, maxLevel); - apf::setScalar(finalSize, v, 0, finalS); - } - m->end(it); -} diff --git a/capstone_clis/capStoneAttachSolution2.cc b/capstone_clis/capStoneAttachSolution2.cc deleted file mode 100644 index 7d73b0592..000000000 --- a/capstone_clis/capStoneAttachSolution2.cc +++ /dev/null @@ -1,1108 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - - -typedef std::vector row; - -std::vector readTable(const char* name); - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col); - -void removeUnusedVerts(apf::Mesh2* m, int offset); - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, - std::vector > &surfToStrandMap); - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize); - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize); - -void writeCre(CapstoneModule& cs, const std::string& fileName); - -struct SortingStruct -{ - apf::Vector3 v; - double wm; - bool operator<(const SortingStruct &other) const - { - return wm < other.wm; - } -}; - -//gradation routines from Proteus - -int gradeSizeModify(apf::Mesh* m, apf::Field* size_iso,double gradingFactor, - double size[2], apf::Adjacent edgAdjVert, - apf::Adjacent vertAdjEdg, - std::queue &markedEdges, - apf::MeshTag* isMarked, - int fieldType, - int vecPos, //which idx of sizeVec to modify - int idxFlag) - -//General function to actually modify sizes -{ - (void)vecPos; - //Determine a switching scheme depending on which vertex needs a modification - int idx1,idx2; - if(idxFlag == 0){ - idx1=0; - idx2=1; - } - else{ - idx1=1; - idx2 = 0; - } - - int marker[3] = {0,1,0}; - double marginVal = 0.01; - int needsParallel=0; - - if(fieldType == apf::SCALAR){ - //apf::Field* size_iso = m->findField("size"); - - if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal)) - { - if(m->isOwned(edgAdjVert[idx1])) - { - size[idx1] = gradingFactor*size[idx2]; - apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]); - m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - //if edge is not already marked - if(!marker[2]){ - m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]); - markedEdges.push(vertAdjEdg[i]); - } - } - } //end isOwned - else - { //Pack information to owning processor - needsParallel=1; - apf::Copies remotes; - m->getRemotes(edgAdjVert[idx1],remotes); - double newSize = gradingFactor*size[idx2]; - int owningPart=m->getOwner(edgAdjVert[idx1]); - PCU_COMM_PACK(owningPart, remotes[owningPart]); - PCU_COMM_PACK(owningPart,newSize); - } - } - - }//end if apf::SCALAR - return needsParallel; -} - -void markEdgesInitial(apf::Mesh* m, apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used to initially determine which edges need to be considered for gradation -{ - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - double size[2]; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::MeshEntity* edge; - apf::MeshIterator* it = m->begin(1); - while((edge=m->iterate(it))){ - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i]=apf::getScalar(size_iso,edgAdjVert[i],0); - } - if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ - //add edge to a queue - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - else{ - m->setIntTag(edge,isMarked,&marker[0]); - } - } - m->end(it); -} - -int serialGradation(apf::Mesh* m,apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used loop over the mesh edge queue for gradation and modify the sizes -{ - double size[2]; - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - apf::MeshEntity* edge; - int needsParallel=0; - - //perform serial gradation while packing necessary info for parallel - while(!markedEdges.empty()){ - edge = markedEdges.front(); - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i] = apf::getScalar(size_iso,edgAdjVert[i],0); - } - - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 0); - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 1); - - m->setIntTag(edge,isMarked,&marker[0]); - markedEdges.pop(); - } - return needsParallel; -} - -int gradeMesh(apf::Mesh* m,apf::Field* size_iso) -//Function to grade isotropic mesh through comparison of edge vertex size ratios -//This implementation accounts for parallel meshes as well -//First do serial gradation. -//If a shared entity has its size modified, then send new size to owning copy. -//After full loop over entities, have owning copy take minimum of all sizes received -//Flag adjacent entities to owning copy. -//Communicate to remote copies that a size was modified, and flag adjacent edges to remote copies for further gradation -{ - //unique to this code - //apf::Field* size_iso = m->findField("size"); - double gradingFactor = 1.2; - // - - apf::MeshEntity* edge; - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - std::queue markedEdges; - apf::MeshTag* isMarked = m->createIntTag("isMarked",1); - - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - apf::MeshIterator* it; - markEdgesInitial(m,size_iso,markedEdges,gradingFactor); - - int needsParallel=1; - while(needsParallel) - { - PCU_Comm_Begin(); - needsParallel = serialGradation(m,size_iso,markedEdges,gradingFactor); - - PCU_Add_Ints(&needsParallel,1); - PCU_Comm_Send(); - - apf::MeshEntity* ent; - double receivedSize; - double currentSize; - double newSize; - - //Need a container to get all entitites that need to be updated on remotes - std::queue updateRemoteVertices; - - apf::Copies remotes; - //owning copies are receiving - while(PCU_Comm_Receive()) - { - PCU_COMM_UNPACK(ent); - PCU_COMM_UNPACK(receivedSize); - - if(!m->isOwned(ent)){ - std::cout<<"THERE WAS AN ERROR"<getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - updateRemoteVertices.push(ent); - } - - PCU_Comm_Begin(); - - while(!updateRemoteVertices.empty()) - { - ent = updateRemoteVertices.front(); - //get remote copies and send updated mesh sizes - m->getRemotes(ent,remotes); - currentSize = apf::getScalar(size_iso,ent,0); - for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter) - { - PCU_COMM_PACK(iter->first, iter->second); - } - updateRemoteVertices.pop(); - } - - PCU_Comm_Send(); - //while remote copies are receiving - while(PCU_Comm_Receive()) - { - //unpack - PCU_COMM_UNPACK(ent); - //PCU_COMM_UNPACK(receivedSize); - assert(!m->isOwned(ent)); - - if(m->isOwned(ent)){ - std::cout<<"Problem occurred\n"; - std::exit(1); - } - - //add adjacent edges into Q - m->getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - } - apf::synchronize(size_iso); - - } //end outer while - - //Cleanup of edge marker field - it = m->begin(1); - while((edge=m->iterate(it))){ - m->removeTag(edge,isMarked); - } - m->end(it); - m->destroyTag(isMarked); - - //apf::synchronize(size_iso); - return needsParallel; -} - -//Eigenvalue routines for volume and surface meshes -// - //function to get volume-based eigenvalues -//apf::Field* getLambdaMax(mesh,hessianField,) -void getLambdaMax(apf::Mesh* mesh,apf::Field* hessianField,apf::Field* lambdaMaxField) -{ - apf::MeshEntity* vert; - apf::MeshIterator* it; - it = mesh->begin(0); - //apf::Field* lambdaMaxField = apf::createLagrangeField(mesh,"lambdaMax",apf::SCALAR,1); - while( (vert = mesh->iterate(it)) ) - { - apf::Matrix3x3 metric; - apf::getMatrix(hessianField, vert, 0, metric); - apf::Vector3 eigenVectors[3]; - double eigenValues[3]; - metric = (metric+apf::transpose(metric))/2.0; - apf::eigen(metric, eigenVectors, eigenValues); - - // Sort the eigenvalues and corresponding vectors - // Larger eigenvalues means a need for a finer mesh - SortingStruct ssa[3]; - for (int i = 0; i < 3; ++i) - { - ssa[i].v = eigenVectors[i]; - ssa[i].wm = std::fabs(eigenValues[i]); - } - std::sort(ssa, ssa + 3); - - assert(ssa[2].wm >= ssa[1].wm); - assert(ssa[1].wm >= ssa[0].wm); - apf::setScalar(lambdaMaxField,vert,0,ssa[2].wm); - } - mesh->end(it); -} - -/* -void getVolMaxPair(apf::Mesh* mesh,apf::Mesh* volMesh,apf::Field* lambdaMaxField,apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshIterator* itVol = volMesh->begin(0); - apf::MeshEntity* vert, *vertVol; - while( (vert = mesh->iterate(it)) ) - { - apf::MeshEntity* vertVol = volMesh->iterate(itVol); - if(lambda_max < apf::getScalar(lambdaMaxField,vertVol,0)) - { - lambda_max = apf::getScalar(lambdaMaxField,vertVol,0); - h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - volMesh->end(itVol); - mesh->end(it); -} -*/ - -void getVolMaxPair(apf::Mesh* mesh,std::vector > surfToStrandMap, apf::Field* lambdaMaxField,apf::Field* lambdaStrandMax, apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert, *vertVol; - int counter = 0; - while( (vert = mesh->iterate(it)) ) - { - //the goal of this loop is to identify lambda max global, find corresponding surface h-max global, and then set the surface lambda field - int strandSize = surfToStrandMap[counter].size(); - for(int i=0; iend(it); -} - -void getSurfMaxPair(apf::Mesh* mesh,apf::Field* lambdaMaxField,apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert; - while( (vert = mesh->iterate(it)) ) - { - if(lambda_max < apf::getScalar(lambdaMaxField,vert,0)) - { - lambda_max = apf::getScalar(lambdaMaxField,vert,0); - h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - mesh->end(it); -} - -void setSizeField(apf::Mesh* mesh, apf::Field* lambdaMaxField,apf::Field* sizeField,double lambda_max,double lambda_cutoff, double h_lambdamax,double h_global,double factor) -{ - apf::MeshIterator* it=mesh->begin(0); - apf::MeshEntity* vert; - int counter3 = 0; - double h_special = -1; - while( (vert = mesh->iterate(it)) ) - { - double h_v = h_special; //set default size as user specified input - double lambda_vert = apf::getScalar(lambdaMaxField,vert,0); - if(lambda_vert > lambda_cutoff) - { - h_v = sqrt(lambda_max/lambda_vert/factor)*h_lambdamax; - } - else - counter3++; - //h_global is user-specified max size - if(h_v > h_global) //maximum value for size field - h_v = h_global; - apf::setScalar(sizeField,vert,0,h_v); - } - mesh->end(it); - -} - -void isotropicIntersect(apf::Mesh* m, std::queue sizeFieldList,apf::Field* finalSizeField,apf::Field* finalChoiceField) -{ - apf::MeshEntity *vert; - apf::MeshIterator *it = m->begin(0); - - apf::Field *field = sizeFieldList.front(); - apf::copyData(finalSizeField,field); - sizeFieldList.pop(); - int choiceIdx = 1; //assumes the initial field was set to choice 0 - while(!sizeFieldList.empty()) - { - field = sizeFieldList.front(); - while(vert = m->iterate(it)) - { - double value1 = apf::getScalar(finalSizeField,vert,0); - double value2 = apf::getScalar(field,vert,0); - double minValue = std::min(value1,value2); - apf::setScalar(finalSizeField,vert,0,minValue); - if(value1 > value2) - { - apf::setScalar(finalChoiceField,vert,0,choiceIdx); - } - } - sizeFieldList.pop(); - choiceIdx++; - } - m->end(it); -} - - -// - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - - double initialTime = PCU_Time(); - - if (argc != 7) { - if(0==PCU_Comm_Self()) - std::cerr << "usage: " << argv[0] - << " \n"; - return EXIT_FAILURE; - } - - - gmi_register_mesh(); - gmi_register_null(); - - const char* creFileName = argv[1]; - const char* dataFileName = argv[2]; - int offset = atoi(argv[3]); - const int strandSize = atoi(argv[4]); - double h_global = atof(argv[5]); - double factor = atof(argv[6]); - - // load capstone mesh - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */ - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(creFileName); - - M_GModel gmodel = cs.load_files(filenames); - - int numbreps = 0; - MG_CALL(g->get_num_breps(numbreps)); - std::cout << "number of b reps is " << numbreps << std::endl; - if(numbreps == 0) - error(HERE, ERR_INVALID_INPUT, "Model is empty"); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - for(std::size_t i = 0; i < mmodels.size(); ++i) - { - M_MModel ammodel = mmodels[i]; - std::size_t numregs = 0; - std::size_t numfaces = 0; - std::size_t numedges = 0; - std::size_t numverts = 0; - MG_API_CALL(m, set_current_model(ammodel)); - MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs)); - MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces)); - MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges)); - MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts)); - std::cout << "num regions is " << numregs << std::endl; - std::cout << "num faces is " << numfaces << std::endl; - std::cout << "num edges is " << numedges << std::endl; - std::cout << "num verts is " << numverts << std::endl; - std::cout << "-----------" << std::endl; - if(numregs > 0) - { - mmodel = ammodel; - break; - } - } - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - - gmi_cap_start(); - gmi_register_cap(); - - printf("---- CapStone Mesh Loaded. \n"); - - // read the data into a vector array for now - printf("\n---- Reading Strand Data. \n"); - std::vector table = readTable(dataFileName); - printf("---- Reading Strand Data: Done. \n"); - - // create the mesh object (this one is CapStone underneath) - printf("\n---- Creating Mesh Wrapper Object. \n"); - apf::Mesh2* mesh = apf::createMesh(m,g); - printf("---- Creating Mesh Wrapper Object: Done. \n"); - - // remove unused verts - if (offset > 0) { - printf("\n---- Removing Extra Verts. \n"); - removeUnusedVerts(mesh, offset); - printf("---- Removing Extra Verts: Done. \n"); - } - - // make the volume mesh (this one is MDS underneath) - printf("\n---- Creating Volume Mesh. \n"); - // create an empty array with the correct sizes to hold the surface 2 strand map - std::vector > surfToStrandMap(mesh->count(0), - std::vector(strandSize, NULL)); - apf::Mesh2* volMesh = createVolumeMesh(mesh, table, strandSize, surfToStrandMap); - printf("---- Creating Volume Mesh: Done. \n"); - - // get all fields and add them to the mesh - printf("\n---- Adding Fields to Volume Mesh. \n"); - apf::Field* rhoField = addScalarField(volMesh, table, "rho", 3, strandSize); - apf::Field* rho_uvwField = addVector3Field(volMesh, table, "rho_uvw", 4, 5, 6, strandSize); - apf::Field* eField = addScalarField(volMesh, table, "e", 7, strandSize); - apf::Field* nuField = addScalarField(volMesh, table, "nu", 8, strandSize); - printf("---- Adding Fields to Volume Mesh: Done. \n"); - - printf("\n---- Printing Volume/Strand Mesh Stats. \n"); - printf("number of mesh verts: %zu\n", volMesh->count(0)); - printf("number of mesh regions(hexes): %zu\n", volMesh->count(3)); - printf("---- Printing Volume/Strand Mesh Stats: Done. \n"); - - double constructionTime = PCU_Time(); - std::cout<<"TIMER: Finished converting capstone mesh to volume mesh "<begin(0); - double minNu = 1.0; - while( (vert = volMesh->iterate(it)) ) - { - double rhoVal = apf::getScalar(rhoField,vert,0); - apf::Vector3 rhoVelocity; - apf::getVector(rho_uvwField,vert,0,rhoVelocity); - double speed = rhoVelocity.getLength()/rhoVal; - apf::setScalar(speedField,vert,0,speed); - if(apf::getScalar(nuField,vert,0)end(it); - apf::Field* gradSpeedField = apf::recoverGradientByVolume(speedField); - apf::Field* hessianSpeedField = apf::recoverGradientByVolume(gradSpeedField); - //End getSpeed - - apf::Field* surfaceSpeedField = apf::createLagrangeField(mesh,"surface_speed",apf::SCALAR,1); - it = mesh->begin(0); - apf::MeshIterator* itVol = volMesh->begin(0); - - while( (vert = mesh->iterate(it)) ) - { - apf::MeshEntity* vertVol = volMesh->iterate(itVol); - double rhoVal = apf::getScalar(rhoField,vertVol,0); - double eVal = apf::getScalar(eField,vertVol,0); - double speedVal = apf::getScalar(speedField,vertVol,0); - apf::setScalar(surfaceRhoField,vert,0,rhoVal); - apf::setScalar(surfaceEField,vert,0,eVal); - apf::setScalar(surfaceSpeedField,vert,0,speedVal); - } - mesh->end(it); - volMesh->end(itVol); - std::cout<<"Finished surface fields from volume\n"; - - double getSurfaceTime = PCU_Time(); - std::cout<<"TIMER: Finished computing speed and transferring fields to surface "<begin(0); - while( (vert = volMesh->iterate(itVol)) ) - { - double rho_e = apf::getScalar(eField,vert,0); - double rho = apf::getScalar(rhoField,vert,0); - double speedVal = apf::getScalar(speedField,vert,0); - double pressure_gamma = rho_e-0.5*rho*speedVal*speedVal; - apf::setScalar(pressureField,vert,0,pressure_gamma); - } - volMesh->end(itVol); - apf::Field* gradPressureField = apf::recoverGradientByVolume(pressureField); - apf::Field* hessianPressureField = apf::recoverGradientByVolume(gradPressureField); - - //get eigenvalues in the volume mesh - getLambdaMax(volMesh,hessianPressureField,lambdaMaxField); - //getLambdaMax(volMesh,hessianEField,lambdaMaxField); - - struct FieldOfInterest{ - FieldOfInterest() - { - lambda_max = 0.0; - h_lambdamax = 0.0; - } - double lambda_max; - double h_lambdamax; - apf::Field* lambdaMaxField; - apf::Field* lambdaStrandMax; - apf::Field* sizeField; - double lambda_cutoff() - { - return lambda_max*1e-10; - } - }; - - - FieldOfInterest eBased; - eBased.lambdaMaxField = lambdaMaxField; - eBased.lambdaStrandMax = apf::createLagrangeField(mesh,"surf_lambda_strandMax",apf::SCALAR,1); - eBased.sizeField = apf::createLagrangeField(mesh,"surface_size_e",apf::SCALAR,1); - - getVolMaxPair(mesh,surfToStrandMap,eBased.lambdaMaxField,eBased.lambdaStrandMax,currentSize,eBased.lambda_max,eBased.h_lambdamax); - - //set size field - - setSizeField(mesh,eBased.lambdaStrandMax,eBased.sizeField,eBased.lambda_max,eBased.lambda_cutoff(),eBased.h_lambdamax,h_global,factor); - - double getPressureTime = PCU_Time(); - std::cout<<"TIMER: Finished pressure size field "<begin(0); - int vIDcounter = 0; //better way to do this would be with a numbering system - //apf::Numbering* nVID = volMesh->findNumbering; //this would be for a volume mesh - - while( (vert = mesh->iterate(it)) ) - { - apf::Vector3 pointVol; - apf::Vector3 pointSurf; - apf::MeshEntity* volVert = surfToStrandMap[vIDcounter][1]; - double speed = apf::getScalar(speedField,volVert,0); - mesh->getPoint(vert,0,pointSurf); - volMesh->getPoint(volVert,0,pointVol); - apf::Vector3 differenceVec = pointVol-pointSurf; - double shearStress = speed/differenceVec.getLength(); - apf::setScalar(surfaceShearField,vert,0,shearStress); - vIDcounter++; - } - mesh->end(it); - - FieldOfInterest shearBased; - shearBased.lambdaMaxField = apf::createLagrangeField(mesh,"surflambdaMax",apf::SCALAR,1); - shearBased.sizeField = apf::createLagrangeField(mesh,"surface_size",apf::SCALAR,1); - - //get surface shear gradient - std::cout<<"Reached gradshearfield\n"; - apf::Field* gradShearField = apf::recoverGradientByVolume(surfaceShearField); - std::cout<<"Got gradshearfield\n"; - - it=mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::Vector3 gradShear; - apf::getVector(gradShearField,vert,0,gradShear); - apf::setScalar(shearBased.lambdaMaxField,vert,0,gradShear.getLength()); - } - mesh->end(it); - - std::cout<<"Got surf lambda max\n"; - getSurfMaxPair(mesh,shearBased.lambdaMaxField,currentSize,shearBased.lambda_max,shearBased.h_lambdamax); - - std::cout<<"Got surf lambda max pair\n"; - //set size field - setSizeField(mesh,shearBased.lambdaMaxField,shearBased.sizeField,shearBased.lambda_max,shearBased.lambda_cutoff(),shearBased.h_lambdamax,h_global,factor); - - std::cout<<"set size field\n"; - - double getShearTime = PCU_Time(); - std::cout<<"TIMER: Finished skin friction size field "< sizeFieldList; - sizeFieldList.push(eBased.sizeField); - sizeFieldList.push(shearBased.sizeField); - apf::Field* finalChoiceField = apf::createLagrangeField(mesh,"finalChoice",apf::SCALAR,1); - it = mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::setScalar(finalChoiceField,vert,0,0); - } - - isotropicIntersect(mesh,sizeFieldList,finalSizeField,finalChoiceField); - - double getIntersectionTime = PCU_Time(); - std::cout<<"TIMER: get mesh intersection "<shouldRunPreZoltan = true; - in->shouldRunMidZoltan = true; - in->shouldRunPostZoltan = true; - in->maximumImbalance = 1.05; - in->maximumIterations = 10; - in->shouldSnap = true; - in->shouldTransferParametric = true; - in->shouldForceAdaptation = true; - ma::adaptVerbose(in); - - double adaptTime = PCU_Time(); - std::cout<<"TIMER: adaptMesh "<countFields();i++) - { - apf::destroyField(mesh->getField(i)); - } - for(int i =0;icountFields();i++) - { - apf::destroyField(volMesh->getField(i)); - } - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} - -std::vector readTable(const char* name) -{ - std::vector table; - std::ifstream file(name); - while(!file.eof()){ - std::string str; - std::getline(file, str); - std::stringstream ss(str); - double value; - row r; - while (ss >> value) - r.push_back(value); - if (r.size() > 0) - table.push_back(r); - } - return table; -} - -static std::vector extractNthStrandData(const std::vector &table, - const int strandSize, const int n, const int col) -{ - std::vector outData; - int i = 0; - while (i*strandSize+(strandSize-n) < int(table.size())) { - outData.push_back(table[i*strandSize+(strandSize - n)][col]); - i++; - } - return outData; -} - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col) -{ - // surface data corresponds to strand 1; - return extractNthStrandData(table, strandSize, 1, col); -} - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::SCALAR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value = t[vid*strandSize+(strandSize - layer)][col]; - apf::setScalar(f, e, 0, value); - count++; - } - m->end(it); - return f; -} - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::VECTOR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value0 = t[vid*strandSize+(strandSize - layer)][col0]; - double value1 = t[vid*strandSize+(strandSize - layer)][col1]; - double value2 = t[vid*strandSize+(strandSize - layer)][col2]; - apf::setVector(f, e, 0, apf::Vector3(value0, value1, value2)); - count++; - } - m->end(it); - return f; -} - - - -void removeUnusedVerts(apf::Mesh2* m, int offset) -{ - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - if (count < offset) { - PCU_ALWAYS_ASSERT(m->countUpward(e) == 0); - m->destroy(e); - } - else - PCU_ALWAYS_ASSERT(m->countUpward(e) != 0); - count++; - } - m->end(it); -} - -static std::vector -readLayerCoordinates(int layer, const std::vector &t, int s) -{ - std::vector xs = extractNthStrandData(t, s, layer, 0); - std::vector ys = extractNthStrandData(t, s, layer, 1); - std::vector zs = extractNthStrandData(t, s, layer, 2); - - std::vector res; - - for (std::size_t i = 0; i < xs.size(); i++) - res.push_back(apf::Vector3(xs[i], ys[i], zs[i])); - - return res; -} - - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, std::vector > &surfToStrandMap) -{ - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "surf_verts_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - // get the original faces - std::vector faces; - it = m->begin(2); - while ( (e = m->iterate(it)) ) - faces.push_back(e); - m->end(it); - - // create and empty mds mesh - apf::Mesh2* vMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 3, false); - apf::Numbering* nLayer = apf::createNumbering(vMesh, "layer_num", apf::getLagrange(1) , 1); - apf::Numbering* nVID = apf::createNumbering(vMesh, "vid_num", apf::getLagrange(1) , 1); - - // process each layer - // layer = 1 corresponds to the surface nodes - // treat it out side of the loop - std::vector coords = readLayerCoordinates(1, t, s); - std::vector lastVs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, 1); - apf::number(nVID, v, 0, 0, i); - lastVs.push_back(v); - } - - PCU_ALWAYS_ASSERT(coords.size() == lastVs.size()); - - // now process the remaining layers - for (int layer = 2; layer < s+1; layer++) { - // read the coords first - std::vector coords = readLayerCoordinates(layer, t, s); - // create the new verts - std::vector vs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, layer); - apf::number(nVID, v, 0, 0, i); - vs.push_back(v); - } - - // create the prisms - for (std::size_t i = 0; i < faces.size(); i++) { - apf::MeshEntity* f = faces[i]; - apf::MeshEntity* dv[3]; - m->getDownward(f, 0, dv); - int id0 = apf::getNumber(n, dv[0], 0, 0); - int id1 = apf::getNumber(n, dv[1], 0, 0); - int id2 = apf::getNumber(n, dv[2], 0, 0); - - apf::MeshEntity* prismVs[6]; - - prismVs[0] = lastVs[id0]; - prismVs[1] = lastVs[id1]; - prismVs[2] = lastVs[id2]; - - prismVs[3] = vs[id0]; - prismVs[4] = vs[id1]; - prismVs[5] = vs[id2]; - - apf::buildElement(vMesh, 0, apf::Mesh::PRISM, prismVs); - } - // set the lastVs to Vs - for (std::size_t i = 0; i < vs.size(); i++) - lastVs[i] = vs[i]; - } - vMesh->acceptChanges(); - apf::deriveMdsModel(vMesh); - apf::verify(vMesh); - - // create the surface to strand map - it = vMesh->begin(0); - while ( (e = vMesh->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0) - 1; - int vid = apf::getNumber(nVID, e, 0, 0); - PCU_ALWAYS_ASSERT((layer >= 0) && (layer < s)); - PCU_ALWAYS_ASSERT((vid >= 0) && (vid < (int)m->count(0))); - surfToStrandMap[vid][layer] = e; - } - vMesh->end(it); - - // return the volume mesh - return vMesh; -} - -void writeCre(CapstoneModule& cs, const std::string& fileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, fileName.c_str(), idmapping); -} diff --git a/capstone_clis/capStoneAttachSolutionStrandOnly.cc b/capstone_clis/capStoneAttachSolutionStrandOnly.cc deleted file mode 100644 index e39994f52..000000000 --- a/capstone_clis/capStoneAttachSolutionStrandOnly.cc +++ /dev/null @@ -1,1161 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - - -typedef std::vector row; - -struct FieldOfInterest{ - FieldOfInterest() - { - lambda_max = 0.0; - h_lambdamax = 0.0; - } - double lambda_max; - double h_lambdamax; - apf::Field* lambdaMaxField; - apf::Field* lambdaStrandMax; - apf::Field* sizeField; - double lambda_cutoff() - { - return lambda_max*1e-10; - } -}; - -std::vector readTable(const char* name); - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col); - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, - std::vector > &surfToStrandMap); - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize); - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize); - -void writeCre(CapstoneModule& cs, const std::string& fileName); - -void writeMdsMesh(apf::Mesh2* m, const char* name, const char* fieldName); - -void computeSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea); - - -void adjustRefinementLevel(apf::Mesh2* m, apf::Field* finalSize, - apf::Field* currentSize, int maxLevel); - -struct SortingStruct -{ - apf::Vector3 v; - double wm; - bool operator<(const SortingStruct &other) const - { - return wm < other.wm; - } -}; - -//gradation routines from Proteus - -int gradeSizeModify(apf::Mesh* m, apf::Field* size_iso,double gradingFactor, - double size[2], apf::Adjacent edgAdjVert, - apf::Adjacent vertAdjEdg, - std::queue &markedEdges, - apf::MeshTag* isMarked, - int fieldType, - int vecPos, //which idx of sizeVec to modify - int idxFlag) - -//General function to actually modify sizes -{ - (void)vecPos; - //Determine a switching scheme depending on which vertex needs a modification - int idx1,idx2; - if(idxFlag == 0){ - idx1=0; - idx2=1; - } - else{ - idx1=1; - idx2 = 0; - } - - int marker[3] = {0,1,0}; - double marginVal = 0.01; - int needsParallel=0; - - if(fieldType == apf::SCALAR){ - //apf::Field* size_iso = m->findField("size"); - - if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal)) - { - if(m->isOwned(edgAdjVert[idx1])) - { - size[idx1] = gradingFactor*size[idx2]; - apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]); - m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - //if edge is not already marked - if(!marker[2]){ - m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]); - markedEdges.push(vertAdjEdg[i]); - } - } - } //end isOwned - else - { //Pack information to owning processor - needsParallel=1; - apf::Copies remotes; - m->getRemotes(edgAdjVert[idx1],remotes); - double newSize = gradingFactor*size[idx2]; - int owningPart=m->getOwner(edgAdjVert[idx1]); - PCU_COMM_PACK(owningPart, remotes[owningPart]); - PCU_COMM_PACK(owningPart,newSize); - } - } - - }//end if apf::SCALAR - return needsParallel; -} - -void markEdgesInitial(apf::Mesh* m, apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used to initially determine which edges need to be considered for gradation -{ - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - double size[2]; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::MeshEntity* edge; - apf::MeshIterator* it = m->begin(1); - while((edge=m->iterate(it))){ - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i]=apf::getScalar(size_iso,edgAdjVert[i],0); - } - if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ - //add edge to a queue - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - else{ - m->setIntTag(edge,isMarked,&marker[0]); - } - } - m->end(it); -} - -int serialGradation(apf::Mesh* m,apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used loop over the mesh edge queue for gradation and modify the sizes -{ - double size[2]; - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - apf::MeshEntity* edge; - int needsParallel=0; - - //perform serial gradation while packing necessary info for parallel - while(!markedEdges.empty()){ - edge = markedEdges.front(); - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i] = apf::getScalar(size_iso,edgAdjVert[i],0); - } - - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 0); - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 1); - - m->setIntTag(edge,isMarked,&marker[0]); - markedEdges.pop(); - } - return needsParallel; -} - -int gradeMesh(apf::Mesh* m,apf::Field* size_iso) -//Function to grade isotropic mesh through comparison of edge vertex size ratios -//This implementation accounts for parallel meshes as well -//First do serial gradation. -//If a shared entity has its size modified, then send new size to owning copy. -//After full loop over entities, have owning copy take minimum of all sizes received -//Flag adjacent entities to owning copy. -//Communicate to remote copies that a size was modified, and flag adjacent edges to remote copies for further gradation -{ - //unique to this code - //apf::Field* size_iso = m->findField("size"); - double gradingFactor = 1.2; - // - - apf::MeshEntity* edge; - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - std::queue markedEdges; - apf::MeshTag* isMarked = m->createIntTag("isMarked",1); - - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - apf::MeshIterator* it; - markEdgesInitial(m,size_iso,markedEdges,gradingFactor); - - int needsParallel=1; - while(needsParallel) - { - PCU_Comm_Begin(); - needsParallel = serialGradation(m,size_iso,markedEdges,gradingFactor); - - PCU_Add_Ints(&needsParallel,1); - PCU_Comm_Send(); - - apf::MeshEntity* ent; - double receivedSize; - double currentSize; - double newSize; - - //Need a container to get all entitites that need to be updated on remotes - std::queue updateRemoteVertices; - - apf::Copies remotes; - //owning copies are receiving - while(PCU_Comm_Receive()) - { - PCU_COMM_UNPACK(ent); - PCU_COMM_UNPACK(receivedSize); - - if(!m->isOwned(ent)){ - std::cout<<"THERE WAS AN ERROR"<getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - updateRemoteVertices.push(ent); - } - - PCU_Comm_Begin(); - - while(!updateRemoteVertices.empty()) - { - ent = updateRemoteVertices.front(); - //get remote copies and send updated mesh sizes - m->getRemotes(ent,remotes); - currentSize = apf::getScalar(size_iso,ent,0); - for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter) - { - PCU_COMM_PACK(iter->first, iter->second); - } - updateRemoteVertices.pop(); - } - - PCU_Comm_Send(); - //while remote copies are receiving - while(PCU_Comm_Receive()) - { - //unpack - PCU_COMM_UNPACK(ent); - //PCU_COMM_UNPACK(receivedSize); - assert(!m->isOwned(ent)); - - if(m->isOwned(ent)){ - std::cout<<"Problem occurred\n"; - std::exit(1); - } - - //add adjacent edges into Q - m->getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - } - apf::synchronize(size_iso); - - } //end outer while - - //Cleanup of edge marker field - it = m->begin(1); - while((edge=m->iterate(it))){ - m->removeTag(edge,isMarked); - } - m->end(it); - m->destroyTag(isMarked); - - //apf::synchronize(size_iso); - return needsParallel; -} - -//Eigenvalue routines for volume and surface meshes -// - //function to get volume-based eigenvalues -//apf::Field* getLambdaMax(mesh,hessianField,) -void getLambdaMax(apf::Mesh* mesh,apf::Field* hessianField,apf::Field* lambdaMaxField) -{ - apf::MeshEntity* vert; - apf::MeshIterator* it; - it = mesh->begin(0); - //apf::Field* lambdaMaxField = apf::createLagrangeField(mesh,"lambdaMax",apf::SCALAR,1); - while( (vert = mesh->iterate(it)) ) - { - apf::Matrix3x3 metric; - apf::getMatrix(hessianField, vert, 0, metric); - apf::Vector3 eigenVectors[3]; - double eigenValues[3]; - metric = (metric+apf::transpose(metric))/2.0; - apf::eigen(metric, eigenVectors, eigenValues); - - // Sort the eigenvalues and corresponding vectors - // Larger eigenvalues means a need for a finer mesh - SortingStruct ssa[3]; - for (int i = 0; i < 3; ++i) - { - ssa[i].v = eigenVectors[i]; - ssa[i].wm = std::fabs(eigenValues[i]); - } - std::sort(ssa, ssa + 3); - - assert(ssa[2].wm >= ssa[1].wm); - assert(ssa[1].wm >= ssa[0].wm); - apf::setScalar(lambdaMaxField,vert,0,ssa[2].wm); - } - mesh->end(it); -} - -void getVolMaxPair(apf::Mesh* mesh,std::vector > surfToStrandMap, apf::Field* lambdaMaxField,apf::Field* lambdaStrandMax, apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert, *vertVol; - int counter = 0; - while( (vert = mesh->iterate(it)) ) - { - //the goal of this loop is to identify lambda max global, find corresponding surface h-max global, and then set the surface lambda field - int strandSize = surfToStrandMap[counter].size(); - for(int i=0; iend(it); -} - -void getSurfMaxPair(apf::Mesh* mesh,apf::Field* lambdaMaxField,apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert; - while( (vert = mesh->iterate(it)) ) - { - if(lambda_max < apf::getScalar(lambdaMaxField,vert,0)) - { - lambda_max = apf::getScalar(lambdaMaxField,vert,0); - h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - mesh->end(it); -} - -void setSizeField(apf::Mesh* mesh, apf::Field* lambdaMaxField,apf::Field* sizeField,double lambda_max,double lambda_cutoff, double h_lambdamax,double h_global,double factor) -{ - apf::MeshIterator* it=mesh->begin(0); - apf::MeshEntity* vert; - int counter3 = 0; - double h_special = -1; - while( (vert = mesh->iterate(it)) ) - { - double h_v = h_special; //set default size as user specified input - double lambda_vert = apf::getScalar(lambdaMaxField,vert,0); - if(lambda_vert > lambda_cutoff) - { - h_v = sqrt(lambda_max/lambda_vert/factor)*h_lambdamax; - } - else - counter3++; - //h_global is user-specified max size - if(h_v > h_global) //maximum value for size field - h_v = h_global; - apf::setScalar(sizeField,vert,0,h_v); - } - mesh->end(it); - -} - -static std::vector decodeBitFields(const char* bitFields) -{ - std::vector res; - res.resize(strlen(bitFields)); - printf("length is %d\n", strlen(bitFields)); - for (int i = 0; i < strlen(bitFields); i++) { - if (bitFields[i] == '0') - res[i] = false; - else if (bitFields[i] == '1') - res[i] = true; - else - PCU_ALWAYS_ASSERT(0); - } - return res; -} - -void isotropicIntersect(apf::Mesh* m, std::queue sizeFieldList, const char* bitFields, apf::Field* finalSizeField,apf::Field* finalChoiceField) -{ - std::vector userInput = decodeBitFields(bitFields); - PCU_ALWAYS_ASSERT(userInput.size() == sizeFieldList.size()); - printf("HERE\n"); - apf::MeshEntity *vert; - apf::MeshIterator *it = m->begin(0); - - apf::Field *field = sizeFieldList.front(); - apf::copyData(finalSizeField,field); - sizeFieldList.pop(); - int choiceIdx = 1; //assumes the initial field was set to choice 0 - while(!sizeFieldList.empty()) - { - field = sizeFieldList.front(); - while( (vert = m->iterate(it)) ) - { - double value1 = apf::getScalar(finalSizeField,vert,0); - double value2 = apf::getScalar(field,vert,0); - if (!userInput[0]) value1 *= 1.e16; - if (!userInput[1]) value2 *= 1.e16; - double minValue = std::min(value1,value2); - apf::setScalar(finalSizeField,vert,0,minValue); - if(value1 > value2) - { - apf::setScalar(finalChoiceField,vert,0,choiceIdx); - } - } - sizeFieldList.pop(); - choiceIdx++; - } - m->end(it); -} - - -// - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - lion_set_verbosity(1); - double initialTime = PCU_Time(); - - if (argc != 7) { - if(0==PCU_Comm_Self()) { - std::cerr << "usage: " << argv[0] - << " \n"; - std::cerr << "*target field(s) is a bit string to select which field(s) are used for error estimation\n"; - std::cerr << "1st bit --- pressure\n"; - std::cerr << "2nd bit --- skin friction\n"; - } - return EXIT_FAILURE; - } - - - gmi_register_mesh(); - gmi_register_null(); - - const char* dataFileName = argv[1]; - const char* bitFields = argv[2]; - const int strandSize = atoi(argv[3]); - double h_global = atof(argv[4]); - double factor = atof(argv[5]); - const int maxLevel = atoi(argv[6]); - - gmi_cap_start(); - gmi_register_cap(); - - // read the data into a vector array for now - printf("\n---- Reading Strand Data. \n"); - std::vector table = readTable(dataFileName); - printf("---- Reading Strand Data: Done. \n"); - - // create the mesh object (this one is CapStone underneath) - printf("\n---- Creating Mesh Wrapper Object. \n"); - apf::Mesh2* mesh = apf::createMesh(m,g); - printf("---- Creating Mesh Wrapper Object: Done. \n"); - - // make the volume mesh (this one is MDS underneath) - printf("\n---- Creating Volume Mesh. \n"); - // create an empty array with the correct sizes to hold the surface 2 strand map - int surfVertCount = - std::vector > surfToStrandMap(surfVertCount, - std::vector(strandSize, NULL)); - apf::Mesh2* volMesh = createVolumeMesh(mesh, table, strandSize, surfToStrandMap); - printf("---- Creating Volume Mesh: Done. \n"); - - // get all fields and add them to the mesh - printf("\n---- Adding Fields to Volume Mesh. \n"); - apf::Field* rhoField = addScalarField(volMesh, table, "rho", 3, strandSize); - apf::Field* rho_uvwField = addVector3Field(volMesh, table, "rho_uvw", 4, 5, 6, strandSize); - apf::Field* eField = addScalarField(volMesh, table, "e", 7, strandSize); - apf::Field* nuField = addScalarField(volMesh, table, "nu", 8, strandSize); - printf("---- Adding Fields to Volume Mesh: Done. \n"); - - printf("\n---- Printing Volume/Strand Mesh Stats. \n"); - printf("number of mesh verts: %zu\n", volMesh->count(0)); - printf("number of mesh regions(hexes): %zu\n", volMesh->count(3)); - printf("---- Printing Volume/Strand Mesh Stats: Done. \n"); - - double constructionTime = PCU_Time(); - std::cout<<"TIMER: Finished converting capstone mesh to volume mesh "<begin(0); - while( (vert = volMesh->iterate(it)) ) - { - double rhoVal = apf::getScalar(rhoField,vert,0); - apf::Vector3 rhoVelocity; - apf::getVector(rho_uvwField,vert,0,rhoVelocity); - double speed = rhoVelocity.getLength()/rhoVal; - apf::setScalar(speedField,vert,0,speed); - } - volMesh->end(it); - //End getSpeed - - apf::Field* surfaceSpeedField = apf::createLagrangeField(mesh,"surface_speed",apf::SCALAR,1); - it = mesh->begin(0); - apf::MeshIterator* itVol = volMesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::MeshEntity* vertVol = volMesh->iterate(itVol); - double rhoVal = apf::getScalar(rhoField,vertVol,0); - double eVal = apf::getScalar(eField,vertVol,0); - double speedVal = apf::getScalar(speedField,vertVol,0); - apf::setScalar(surfaceRhoField,vert,0,rhoVal); - apf::setScalar(surfaceEField,vert,0,eVal); - apf::setScalar(surfaceSpeedField,vert,0,speedVal); - } - mesh->end(it); - volMesh->end(itVol); - std::cout<<"Finished surface fields from volume\n"; - - double getSurfaceTime = PCU_Time(); - std::cout<<"TIMER: Finished computing speed and transferring fields to surface "<begin(0); - while( (vert = volMesh->iterate(itVol)) ) - { - double rho_e = apf::getScalar(eField,vert,0); - double rho = apf::getScalar(rhoField,vert,0); - double speedVal = apf::getScalar(speedField,vert,0); - double pressure_gamma = rho_e-0.5*rho*speedVal*speedVal; - apf::setScalar(pressureField,vert,0,pressure_gamma); - } - volMesh->end(itVol); - apf::Field* gradPressureField = apf::recoverGradientByVolume(pressureField); - apf::Field* hessianPressureField = apf::recoverGradientByVolume(gradPressureField); - - //get eigenvalues in the volume mesh - getLambdaMax(volMesh,hessianPressureField,lambdaMaxField); - - - FieldOfInterest eBased; - eBased.lambdaMaxField = lambdaMaxField; - eBased.lambdaStrandMax = apf::createLagrangeField(mesh,"surf_lambda_strandMax",apf::SCALAR,1); - eBased.sizeField = apf::createLagrangeField(mesh,"surface_size_e",apf::SCALAR,1); - - getVolMaxPair(mesh,surfToStrandMap,eBased.lambdaMaxField,eBased.lambdaStrandMax,currentSize,eBased.lambda_max,eBased.h_lambdamax); - - //set size field - - setSizeField(mesh,eBased.lambdaStrandMax,eBased.sizeField,eBased.lambda_max,eBased.lambda_cutoff(),eBased.h_lambdamax,h_global,factor); - - double getPressureTime = PCU_Time(); - std::cout<<"TIMER: Finished pressure size field "<begin(0); - int vIDcounter = 0; //better way to do this would be with a numbering system - //apf::Numbering* nVID = volMesh->findNumbering; //this would be for a volume mesh - - while( (vert = mesh->iterate(it)) ) - { - apf::Vector3 pointVol; - apf::Vector3 pointSurf; - apf::MeshEntity* volVert = surfToStrandMap[vIDcounter][1]; - double speed = apf::getScalar(speedField,volVert,0); - mesh->getPoint(vert,0,pointSurf); - volMesh->getPoint(volVert,0,pointVol); - apf::Vector3 differenceVec = pointVol-pointSurf; - double shearStress = speed/differenceVec.getLength(); - apf::setScalar(surfaceShearField,vert,0,shearStress); - vIDcounter++; - } - mesh->end(it); - - FieldOfInterest shearBased; - shearBased.lambdaMaxField = apf::createLagrangeField(mesh,"surflambdaMax",apf::SCALAR,1); - shearBased.sizeField = apf::createLagrangeField(mesh,"surface_size",apf::SCALAR,1); - - //get surface shear gradient - std::cout<<"Reached gradshearfield\n"; - apf::Field* gradShearField = apf::recoverGradientByVolume(surfaceShearField); - std::cout<<"Got gradshearfield\n"; - - it=mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::Vector3 gradShear; - apf::getVector(gradShearField,vert,0,gradShear); - apf::setScalar(shearBased.lambdaMaxField,vert,0,gradShear.getLength()); - } - mesh->end(it); - - std::cout<<"Got surf lambda max\n"; - getSurfMaxPair(mesh,shearBased.lambdaMaxField,currentSize,shearBased.lambda_max,shearBased.h_lambdamax); - - std::cout<<"Got surf lambda max pair\n"; - //set size field - setSizeField(mesh,shearBased.lambdaMaxField,shearBased.sizeField,shearBased.lambda_max,shearBased.lambda_cutoff(),shearBased.h_lambdamax,h_global,factor); - - std::cout<<"set size field\n"; - - double getShearTime = PCU_Time(); - std::cout<<"TIMER: Finished skin friction size field "< sizeFieldList; - sizeFieldList.push(eBased.sizeField); - sizeFieldList.push(shearBased.sizeField); - apf::Field* finalChoiceField = apf::createLagrangeField(mesh,"finalChoice",apf::SCALAR,1); - it = mesh->begin(0); - while( (vert = mesh->iterate(it)) ) - { - apf::setScalar(finalChoiceField,vert,0,0); - } - - isotropicIntersect(mesh,sizeFieldList,bitFields,finalSizeField,finalChoiceField); - - double getIntersectionTime = PCU_Time(); - std::cout<<"TIMER: get mesh intersection "<shouldRunPreZoltan = true; - in->shouldRunMidZoltan = true; - in->shouldRunPostZoltan = true; - in->maximumImbalance = 1.05; - in->maximumIterations = log2(factor) + 2; - in->shouldSnap = true; - in->shouldTransferParametric = true; - /* in->shouldTransferToClosestPoint = true; */ - in->shouldForceAdaptation = true; - /* in->debugFolder = "debug"; */ - ma::adaptVerbose(in, false); - - double adaptTime = PCU_Time(); - std::cout<<"TIMER: adaptMesh "< binCount; - std::vector binArea; - computeSizeDistribution(mesh, factor, binCount, binArea); - - for (std::size_t i = 0; i < binCount.size(); i++) { - printf("bin %lu's count/area is %d/%g \n", i, binCount[i], binArea[i]); - } - - - // write the adapted mesh to vtk - apf::writeVtkFiles("adaptedMesh", mesh); - // write the adapted mesh to smb - // Note: This is for statistic measurements using - // measureIsoStats.cc - // The last argument is the name of the field to be used - // in measureIsoStats.cc - writeMdsMesh(mesh, "adaptedMesh.smb", "adapt_size"); - // write the adapted mesh to cre - writeCre(cs, "adaptedMesh.cre"); - - // clean up and exit calls - //destroy all fields... - for(int i =0;icountFields();i++) - { - apf::destroyField(mesh->getField(i)); - } - for(int i =0;icountFields();i++) - { - apf::destroyField(volMesh->getField(i)); - } - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} - -std::vector readTable(const char* name) -{ - std::vector table; - std::ifstream file(name); - while(!file.eof()){ - std::string str; - std::getline(file, str); - std::stringstream ss(str); - double value; - row r; - while (ss >> value) - r.push_back(value); - if (r.size() > 0) - table.push_back(r); - } - return table; -} - -static std::vector extractNthStrandData(const std::vector &table, - const int strandSize, const int n, const int col) -{ - std::vector outData; - int i = 0; - while (i*strandSize+(strandSize-n) < int(table.size())) { - outData.push_back(table[i*strandSize+(strandSize - n)][col]); - i++; - } - return outData; -} - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col) -{ - // surface data corresponds to strand 1; - return extractNthStrandData(table, strandSize, 1, col); -} - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::SCALAR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value = t[vid*strandSize+(strandSize - layer)][col]; - apf::setScalar(f, e, 0, value); - count++; - } - m->end(it); - return f; -} - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::VECTOR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value0 = t[vid*strandSize+(strandSize - layer)][col0]; - double value1 = t[vid*strandSize+(strandSize - layer)][col1]; - double value2 = t[vid*strandSize+(strandSize - layer)][col2]; - apf::setVector(f, e, 0, apf::Vector3(value0, value1, value2)); - count++; - } - m->end(it); - return f; -} - - - -static std::vector -readLayerCoordinates(int layer, const std::vector &t, int s) -{ - std::vector xs = extractNthStrandData(t, s, layer, 0); - std::vector ys = extractNthStrandData(t, s, layer, 1); - std::vector zs = extractNthStrandData(t, s, layer, 2); - - std::vector res; - - for (std::size_t i = 0; i < xs.size(); i++) - res.push_back(apf::Vector3(xs[i], ys[i], zs[i])); - - return res; -} - - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, std::vector > &surfToStrandMap) -{ - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "surf_verts_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - // get the original faces - std::vector faces; - it = m->begin(2); - while ( (e = m->iterate(it)) ) - faces.push_back(e); - m->end(it); - - // create and empty mds mesh - apf::Mesh2* vMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 3, false); - apf::Numbering* nLayer = apf::createNumbering(vMesh, "layer_num", apf::getLagrange(1) , 1); - apf::Numbering* nVID = apf::createNumbering(vMesh, "vid_num", apf::getLagrange(1) , 1); - - // process each layer - // layer = 1 corresponds to the surface nodes - // treat it out side of the loop - std::vector coords = readLayerCoordinates(1, t, s); - std::vector lastVs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, 1); - apf::number(nVID, v, 0, 0, i); - lastVs.push_back(v); - } - - PCU_ALWAYS_ASSERT(coords.size() == lastVs.size()); - - // now process the remaining layers - for (int layer = 2; layer < s+1; layer++) { - // read the coords first - std::vector coords = readLayerCoordinates(layer, t, s); - // create the new verts - std::vector vs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, layer); - apf::number(nVID, v, 0, 0, i); - vs.push_back(v); - } - - // create the prisms - for (std::size_t i = 0; i < faces.size(); i++) { - apf::MeshEntity* f = faces[i]; - apf::MeshEntity* dv[3]; - m->getDownward(f, 0, dv); - int id0 = apf::getNumber(n, dv[0], 0, 0); - int id1 = apf::getNumber(n, dv[1], 0, 0); - int id2 = apf::getNumber(n, dv[2], 0, 0); - - apf::MeshEntity* prismVs[6]; - - prismVs[0] = lastVs[id0]; - prismVs[1] = lastVs[id1]; - prismVs[2] = lastVs[id2]; - - prismVs[3] = vs[id0]; - prismVs[4] = vs[id1]; - prismVs[5] = vs[id2]; - - apf::buildElement(vMesh, 0, apf::Mesh::PRISM, prismVs); - } - // set the lastVs to Vs - for (std::size_t i = 0; i < vs.size(); i++) - lastVs[i] = vs[i]; - } - vMesh->acceptChanges(); - apf::deriveMdsModel(vMesh); - apf::verify(vMesh); - - // create the surface to strand map - it = vMesh->begin(0); - while ( (e = vMesh->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0) - 1; - int vid = apf::getNumber(nVID, e, 0, 0); - PCU_ALWAYS_ASSERT((layer >= 0) && (layer < s)); - PCU_ALWAYS_ASSERT((vid >= 0) && (vid < (int)m->count(0))); - surfToStrandMap[vid][layer] = e; - } - vMesh->end(it); - - // return the volume mesh - return vMesh; -} - -void writeCre(CapstoneModule& cs, const std::string& fileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, fileName.c_str(), idmapping); -} - -void writeMdsMesh(apf::Mesh2* m, const char* name, const char* fieldName) -{ - // this would be the smb mesh with a field associated with it - apf::Mesh2* outMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 2, false); - - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "v_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - int numVerts = m->count(0); - apf::MeshEntity* vMap[numVerts]; - - apf::MeshEntity* v; - it = m->begin(0); - while( (v = m->iterate(it)) ) { - apf::Vector3 p; - m->getPoint(v, 0, p); - apf::MeshEntity* nv = outMesh->createVertex(0, p, apf::Vector3(0, 0, 0)); - int id = apf::getNumber(n, v, 0, 0); - vMap[id] = nv; - } - m->end(it); - - it = m->begin(2); - while( (e = m->iterate(it)) ) { - apf::MeshEntity* vs[3]; - m->getDownward(e, 0, vs); - apf::MeshEntity* triVs[3]; - triVs[0] = vMap[apf::getNumber(n, vs[0], 0, 0)]; - triVs[1] = vMap[apf::getNumber(n, vs[1], 0, 0)]; - triVs[2] = vMap[apf::getNumber(n, vs[2], 0, 0)]; - apf::buildElement(outMesh, 0, apf::Mesh::TRIANGLE, triVs); - } - m->end(it); - - outMesh->acceptChanges(); - apf::deriveMdsModel(outMesh); - outMesh->verify(); - apf::verify(outMesh); - - apf::destroyNumbering(n); - - // get the field on input mesh - apf::Field* inputField = m->findField(fieldName); - PCU_ALWAYS_ASSERT(inputField); - - // create the corresponding field for the out mesh - apf::Field* outputField = apf::createFieldOn(outMesh, "sizes", apf::SCALAR); - - it = m->begin(0); - apf::MeshIterator* outIt = outMesh->begin(0); - while( (v = m->iterate(it)) ) { - double value = apf::getScalar(inputField, v, 0); - apf::MeshEntity* outV = outMesh->iterate(outIt); - apf::setScalar(outputField, outV, 0, value); - } - m->end(it); - outMesh->end(outIt); - outMesh->writeNative(name); - - // clean up outMesh - outMesh->destroyNative(); - apf::destroyMesh(outMesh); -} - - -void computeSizeDistribution(apf::Mesh2* m, int factor, - std::vector& binCount, std::vector& binArea) -{ - // find the min length edge - apf::MeshEntity* e; - apf::MeshIterator* it; - - double minLength = 1.0e12; - it = m->begin(1); - while ( (e = m->iterate(it)) ) { - double l = apf::measure(m, e); - if (l < minLength) - minLength = l; - } - m->end(it); - - double minArea = minLength * minLength / 2.; - - // initialize the arrays - for (int i = 1; i <= factor; i*=2) { - binCount.push_back(0); - binArea.push_back(0.0); - } - - - it = m->begin(2); - while ( (e = m->iterate(it)) ) { - double a = apf::measure(m, e); - int count = 0; - for (int i = 1; i <= factor; i*=2) { - if (a >= i * minArea && a < 4 * i * minArea) { - binCount[count] += 1; - binArea[count] += a; - break; - } - count++; - } - } - m->end(it); -} - -void adjustRefinementLevel(apf::Mesh2* m, apf::Field* finalSize, - apf::Field* currentSize, int maxLevel) -{ - apf::MeshEntity* v; - apf::MeshIterator* it = m->begin(0); - - while( (v = m->iterate(it)) ) { - double currnetS = apf::getScalar(currentSize, v, 0); - double finalS = apf::getScalar(finalSize, v, 0); - if (currnetS < finalS) continue; - if (log2(currnetS/finalS) > maxLevel) - finalS = currnetS / pow(2, maxLevel); - apf::setScalar(finalSize, v, 0, finalS); - } - m->end(it); -} diff --git a/capstone_clis/capStoneAttachSolutionUni.cc b/capstone_clis/capStoneAttachSolutionUni.cc deleted file mode 100644 index 504cc9e03..000000000 --- a/capstone_clis/capStoneAttachSolutionUni.cc +++ /dev/null @@ -1,951 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - - -typedef std::vector row; - -std::vector readTable(const char* name); - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col); - -void removeUnusedVerts(apf::Mesh2* m, int offset); - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, - std::vector > &surfToStrandMap); - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize); - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize); - -void writeCre(CapstoneModule& cs, const std::string& fileName); - -void writeMdsMesh(apf::Mesh2* m, const char* name, const char* fieldName); - -struct SortingStruct -{ - apf::Vector3 v; - double wm; - bool operator<(const SortingStruct &other) const - { - return wm < other.wm; - } -}; - -//gradation routines from Proteus - -int gradeSizeModify(apf::Mesh* m, apf::Field* size_iso,double gradingFactor, - double size[2], apf::Adjacent edgAdjVert, - apf::Adjacent vertAdjEdg, - std::queue &markedEdges, - apf::MeshTag* isMarked, - int fieldType, - int vecPos, //which idx of sizeVec to modify - int idxFlag) - -//General function to actually modify sizes -{ - (void)vecPos; - //Determine a switching scheme depending on which vertex needs a modification - int idx1,idx2; - if(idxFlag == 0){ - idx1=0; - idx2=1; - } - else{ - idx1=1; - idx2 = 0; - } - - int marker[3] = {0,1,0}; - double marginVal = 0.01; - int needsParallel=0; - - if(fieldType == apf::SCALAR){ - //apf::Field* size_iso = m->findField("size"); - - if(size[idx1]>(gradingFactor*size[idx2])*(1+marginVal)) - { - if(m->isOwned(edgAdjVert[idx1])) - { - size[idx1] = gradingFactor*size[idx2]; - apf::setScalar(size_iso,edgAdjVert[idx1],0,size[idx1]); - m->getAdjacent(edgAdjVert[idx1], 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - //if edge is not already marked - if(!marker[2]){ - m->setIntTag(vertAdjEdg[i],isMarked,&marker[1]); - markedEdges.push(vertAdjEdg[i]); - } - } - } //end isOwned - else - { //Pack information to owning processor - needsParallel=1; - apf::Copies remotes; - m->getRemotes(edgAdjVert[idx1],remotes); - double newSize = gradingFactor*size[idx2]; - int owningPart=m->getOwner(edgAdjVert[idx1]); - PCU_COMM_PACK(owningPart, remotes[owningPart]); - PCU_COMM_PACK(owningPart,newSize); - } - } - - }//end if apf::SCALAR - return needsParallel; -} - -void markEdgesInitial(apf::Mesh* m, apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used to initially determine which edges need to be considered for gradation -{ - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - double size[2]; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::MeshEntity* edge; - apf::MeshIterator* it = m->begin(1); - while((edge=m->iterate(it))){ - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i]=apf::getScalar(size_iso,edgAdjVert[i],0); - } - if( (size[0] > gradingFactor*size[1]) || (size[1] > gradingFactor*size[0]) ){ - //add edge to a queue - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - else{ - m->setIntTag(edge,isMarked,&marker[0]); - } - } - m->end(it); -} - -int serialGradation(apf::Mesh* m,apf::Field* size_iso, std::queue &markedEdges,double gradingFactor) -//Function used loop over the mesh edge queue for gradation and modify the sizes -{ - double size[2]; - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - apf::MeshTag* isMarked = m->findTag("isMarked"); -// apf::Field* size_iso = m->findField("size"); - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - apf::MeshEntity* edge; - int needsParallel=0; - - //perform serial gradation while packing necessary info for parallel - while(!markedEdges.empty()){ - edge = markedEdges.front(); - m->getAdjacent(edge, 0, edgAdjVert); - for (std::size_t i=0; i < edgAdjVert.getSize(); ++i){ - size[i] = apf::getScalar(size_iso,edgAdjVert[i],0); - } - - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 0); - needsParallel+=gradeSizeModify(m,size_iso, gradingFactor, size, edgAdjVert, - vertAdjEdg, markedEdges, isMarked, apf::SCALAR,0, 1); - - m->setIntTag(edge,isMarked,&marker[0]); - markedEdges.pop(); - } - return needsParallel; -} - -int gradeMesh(apf::Mesh* m,apf::Field* size_iso) -//Function to grade isotropic mesh through comparison of edge vertex size ratios -//This implementation accounts for parallel meshes as well -//First do serial gradation. -//If a shared entity has its size modified, then send new size to owning copy. -//After full loop over entities, have owning copy take minimum of all sizes received -//Flag adjacent entities to owning copy. -//Communicate to remote copies that a size was modified, and flag adjacent edges to remote copies for further gradation -{ - //unique to this code - //apf::Field* size_iso = m->findField("size"); - double gradingFactor = 1.2; - // - - apf::MeshEntity* edge; - apf::Adjacent edgAdjVert; - apf::Adjacent vertAdjEdg; - std::queue markedEdges; - apf::MeshTag* isMarked = m->createIntTag("isMarked",1); - - //marker structure for 0) not marked 1) marked 2)storage - int marker[3] = {0,1,0}; - - apf::MeshIterator* it; - markEdgesInitial(m,size_iso,markedEdges,gradingFactor); - - int needsParallel=1; - while(needsParallel) - { - PCU_Comm_Begin(); - needsParallel = serialGradation(m,size_iso,markedEdges,gradingFactor); - - PCU_Add_Ints(&needsParallel,1); - PCU_Comm_Send(); - - apf::MeshEntity* ent; - double receivedSize; - double currentSize; - double newSize; - - //Need a container to get all entitites that need to be updated on remotes - std::queue updateRemoteVertices; - - apf::Copies remotes; - //owning copies are receiving - while(PCU_Comm_Receive()) - { - PCU_COMM_UNPACK(ent); - PCU_COMM_UNPACK(receivedSize); - - if(!m->isOwned(ent)){ - std::cout<<"THERE WAS AN ERROR"<getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - updateRemoteVertices.push(ent); - } - - PCU_Comm_Begin(); - - while(!updateRemoteVertices.empty()) - { - ent = updateRemoteVertices.front(); - //get remote copies and send updated mesh sizes - m->getRemotes(ent,remotes); - currentSize = apf::getScalar(size_iso,ent,0); - for(apf::Copies::iterator iter=remotes.begin(); iter!=remotes.end();++iter) - { - PCU_COMM_PACK(iter->first, iter->second); - } - updateRemoteVertices.pop(); - } - - PCU_Comm_Send(); - //while remote copies are receiving - while(PCU_Comm_Receive()) - { - //unpack - PCU_COMM_UNPACK(ent); - //PCU_COMM_UNPACK(receivedSize); - assert(!m->isOwned(ent)); - - if(m->isOwned(ent)){ - std::cout<<"Problem occurred\n"; - std::exit(1); - } - - //add adjacent edges into Q - m->getAdjacent(ent, 1, vertAdjEdg); - for (std::size_t i=0; igetIntTag(vertAdjEdg[i],isMarked,&marker[2]); - if(!marker[2]) - { - markedEdges.push(edge); - //tag edge to indicate that it is part of queue - m->setIntTag(edge,isMarked,&marker[1]); - } - } - } - apf::synchronize(size_iso); - - } //end outer while - - //Cleanup of edge marker field - it = m->begin(1); - while((edge=m->iterate(it))){ - m->removeTag(edge,isMarked); - } - m->end(it); - m->destroyTag(isMarked); - - //apf::synchronize(size_iso); - return needsParallel; -} - -//Eigenvalue routines for volume and surface meshes -// - //function to get volume-based eigenvalues -//apf::Field* getLambdaMax(mesh,hessianField,) -void getLambdaMax(apf::Mesh* mesh,apf::Field* hessianField,apf::Field* lambdaMaxField) -{ - apf::MeshEntity* vert; - apf::MeshIterator* it; - it = mesh->begin(0); - //apf::Field* lambdaMaxField = apf::createLagrangeField(mesh,"lambdaMax",apf::SCALAR,1); - while( (vert = mesh->iterate(it)) ) - { - apf::Matrix3x3 metric; - apf::getMatrix(hessianField, vert, 0, metric); - apf::Vector3 eigenVectors[3]; - double eigenValues[3]; - metric = (metric+apf::transpose(metric))/2.0; - apf::eigen(metric, eigenVectors, eigenValues); - - // Sort the eigenvalues and corresponding vectors - // Larger eigenvalues means a need for a finer mesh - SortingStruct ssa[3]; - for (int i = 0; i < 3; ++i) - { - ssa[i].v = eigenVectors[i]; - ssa[i].wm = std::fabs(eigenValues[i]); - } - std::sort(ssa, ssa + 3); - - assert(ssa[2].wm >= ssa[1].wm); - assert(ssa[1].wm >= ssa[0].wm); - apf::setScalar(lambdaMaxField,vert,0,ssa[2].wm); - } - mesh->end(it); -} - -/* -void getVolMaxPair(apf::Mesh* mesh,apf::Mesh* volMesh,apf::Field* lambdaMaxField,apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshIterator* itVol = volMesh->begin(0); - apf::MeshEntity* vert, *vertVol; - while( (vert = mesh->iterate(it)) ) - { - apf::MeshEntity* vertVol = volMesh->iterate(itVol); - if(lambda_max < apf::getScalar(lambdaMaxField,vertVol,0)) - { - lambda_max = apf::getScalar(lambdaMaxField,vertVol,0); - h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - volMesh->end(itVol); - mesh->end(it); -} -*/ - -void getVolMaxPair(apf::Mesh* mesh,std::vector > surfToStrandMap, apf::Field* lambdaMaxField,apf::Field* lambdaStrandMax, apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert, *vertVol; - int counter = 0; - while( (vert = mesh->iterate(it)) ) - { - //the goal of this loop is to identify lambda max global, find corresponding surface h-max global, and then set the surface lambda field - int strandSize = surfToStrandMap[counter].size(); - for(int i=0; iend(it); -} - -void getSurfMaxPair(apf::Mesh* mesh,apf::Field* lambdaMaxField,apf::Field* currentSize,double &lambda_max,double &h_lambdamax) -{ - apf::MeshIterator* it = mesh->begin(0); - apf::MeshEntity* vert; - while( (vert = mesh->iterate(it)) ) - { - if(lambda_max < apf::getScalar(lambdaMaxField,vert,0)) - { - lambda_max = apf::getScalar(lambdaMaxField,vert,0); - h_lambdamax = apf::getScalar(currentSize,vert,0); - } - } - mesh->end(it); -} - -void setSizeField(apf::Mesh* mesh, apf::Field* lambdaMaxField,apf::Field* sizeField,double lambda_max,double lambda_cutoff, double h_lambdamax,double h_global,double factor) -{ - apf::MeshIterator* it=mesh->begin(0); - apf::MeshEntity* vert; - int counter3 = 0; - double h_special = -1; - while( (vert = mesh->iterate(it)) ) - { - double h_v = h_special; //set default size as user specified input - double lambda_vert = apf::getScalar(lambdaMaxField,vert,0); - if(lambda_vert > lambda_cutoff) - { - h_v = sqrt(lambda_max/lambda_vert/factor)*h_lambdamax; - } - else - counter3++; - //h_global is user-specified max size - if(h_v > h_global) //maximum value for size field - h_v = h_global; - apf::setScalar(sizeField,vert,0,h_v); - } - mesh->end(it); - -} - -void isotropicIntersect(apf::Mesh* m, std::queue sizeFieldList,apf::Field* finalSizeField,apf::Field* finalChoiceField) -{ - apf::MeshEntity *vert; - apf::MeshIterator *it = m->begin(0); - - apf::Field *field = sizeFieldList.front(); - apf::copyData(finalSizeField,field); - sizeFieldList.pop(); - int choiceIdx = 1; //assumes the initial field was set to choice 0 - while(!sizeFieldList.empty()) - { - field = sizeFieldList.front(); - while( (vert = m->iterate(it)) ) - { - double value1 = apf::getScalar(finalSizeField,vert,0); - double value2 = apf::getScalar(field,vert,0); - double minValue = std::min(value1,value2); - apf::setScalar(finalSizeField,vert,0,minValue); - if(value1 > value2) - { - apf::setScalar(finalChoiceField,vert,0,choiceIdx); - } - } - sizeFieldList.pop(); - choiceIdx++; - } - m->end(it); -} - - -// - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - lion_set_verbosity(1); - double initialTime = PCU_Time(); - - if (argc != 7) { - if(0==PCU_Comm_Self()) - std::cerr << "usage: " << argv[0] - << " \n"; - return EXIT_FAILURE; - } - - - gmi_register_mesh(); - gmi_register_null(); - - const char* creFileName = argv[1]; - const char* dataFileName = argv[2]; - int offset = atoi(argv[3]); - const int strandSize = atoi(argv[4]); - double h_global = atof(argv[5]); - double factor = atof(argv[6]); - - // load capstone mesh - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - /* const std::string gdbName("Geometry Database : Create");// Switch Create with SMLIB for CAD */ - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("the_module", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(creFileName); - - M_GModel gmodel = cs.load_files(filenames); - - int numbreps = 0; - MG_CALL(g->get_num_breps(numbreps)); - std::cout << "number of b reps is " << numbreps << std::endl; - if(numbreps == 0) - error(HERE, ERR_INVALID_INPUT, "Model is empty"); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - for(std::size_t i = 0; i < mmodels.size(); ++i) - { - M_MModel ammodel = mmodels[i]; - std::size_t numregs = 0; - std::size_t numfaces = 0; - std::size_t numedges = 0; - std::size_t numverts = 0; - MG_API_CALL(m, set_current_model(ammodel)); - MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs)); - MG_API_CALL(m, get_num_topos(TOPO_FACE, numfaces)); - MG_API_CALL(m, get_num_topos(TOPO_EDGE, numedges)); - MG_API_CALL(m, get_num_topos(TOPO_VERTEX, numverts)); - std::cout << "num regions is " << numregs << std::endl; - std::cout << "num faces is " << numfaces << std::endl; - std::cout << "num edges is " << numedges << std::endl; - std::cout << "num verts is " << numverts << std::endl; - std::cout << "-----------" << std::endl; - if(numregs > 0) - { - mmodel = ammodel; - break; - } - } - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - - gmi_cap_start(); - gmi_register_cap(); - - printf("---- CapStone Mesh Loaded. \n"); - - apf::Mesh2* mesh = apf::createMesh(m,g); - - //adapt the mesh - ma::Input* in; - in = ma::configureUniformRefine(mesh, 1); - ma::validateInput(in); - in->shouldRunPreZoltan = true; - in->shouldRunMidZoltan = true; - in->shouldRunPostZoltan = true; - in->maximumImbalance = 1.05; - in->maximumIterations = log2(factor) + 2; - in->shouldSnap = true; - in->shouldTransferParametric = true; - /* in->shouldTransferToClosestPoint = true; */ - in->shouldForceAdaptation = true; - in->debugFolder = "debug"; - ma::adaptVerbose(in, true); - - double adaptTime = PCU_Time(); - - - // write the adapted mesh to vtk - apf::writeVtkFiles("adaptedMesh", mesh); - // write the adapted mesh to smb - // Note: This is for statistic measurements using - // measureIsoStats.cc - // The last argument is the name of the field to be used - // in measureIsoStats.cc - writeMdsMesh(mesh, "adaptedMesh.smb", "adapt_size"); - // write the adapted mesh to cre - writeCre(cs, "adaptedMesh.cre"); - - - // clean up and exit calls - //destroy all fields... - for(int i =0;icountFields();i++) - { - apf::destroyField(mesh->getField(i)); - } - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} - -std::vector readTable(const char* name) -{ - std::vector table; - std::ifstream file(name); - while(!file.eof()){ - std::string str; - std::getline(file, str); - std::stringstream ss(str); - double value; - row r; - while (ss >> value) - r.push_back(value); - if (r.size() > 0) - table.push_back(r); - } - return table; -} - -static std::vector extractNthStrandData(const std::vector &table, - const int strandSize, const int n, const int col) -{ - std::vector outData; - int i = 0; - while (i*strandSize+(strandSize-n) < int(table.size())) { - outData.push_back(table[i*strandSize+(strandSize - n)][col]); - i++; - } - return outData; -} - -std::vector extractSurfaceData(const std::vector &table, - const int strandSize, const int col) -{ - // surface data corresponds to strand 1; - return extractNthStrandData(table, strandSize, 1, col); -} - -apf::Field* addScalarField(apf::Mesh2* m, const std::vector t, const char* name, int col, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::SCALAR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value = t[vid*strandSize+(strandSize - layer)][col]; - apf::setScalar(f, e, 0, value); - count++; - } - m->end(it); - return f; -} - -apf::Field* addVector3Field(apf::Mesh2* m, const std::vector t, const char* name, - int col0, int col1, int col2, int strandSize) -{ - apf::Numbering* nLayer = m->findNumbering("layer_num"); - apf::Numbering* nVid = m->findNumbering("vid_num"); - PCU_ALWAYS_ASSERT(nLayer); - PCU_ALWAYS_ASSERT(nVid); - - apf::Field* f = apf::createFieldOn(m, name, apf::VECTOR); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0); - int vid = apf::getNumber(nVid, e, 0, 0); - double value0 = t[vid*strandSize+(strandSize - layer)][col0]; - double value1 = t[vid*strandSize+(strandSize - layer)][col1]; - double value2 = t[vid*strandSize+(strandSize - layer)][col2]; - apf::setVector(f, e, 0, apf::Vector3(value0, value1, value2)); - count++; - } - m->end(it); - return f; -} - - - -void removeUnusedVerts(apf::Mesh2* m, int offset) -{ - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - if (count < offset) { - PCU_ALWAYS_ASSERT(m->countUpward(e) == 0); - m->destroy(e); - } - else - PCU_ALWAYS_ASSERT(m->countUpward(e) != 0); - count++; - } - m->end(it); -} - -static std::vector -readLayerCoordinates(int layer, const std::vector &t, int s) -{ - std::vector xs = extractNthStrandData(t, s, layer, 0); - std::vector ys = extractNthStrandData(t, s, layer, 1); - std::vector zs = extractNthStrandData(t, s, layer, 2); - - std::vector res; - - for (std::size_t i = 0; i < xs.size(); i++) - res.push_back(apf::Vector3(xs[i], ys[i], zs[i])); - - return res; -} - - -apf::Mesh2* createVolumeMesh(apf::Mesh2* m, const std::vector &t, int s, std::vector > &surfToStrandMap) -{ - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "surf_verts_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - // get the original faces - std::vector faces; - it = m->begin(2); - while ( (e = m->iterate(it)) ) - faces.push_back(e); - m->end(it); - - // create and empty mds mesh - apf::Mesh2* vMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 3, false); - apf::Numbering* nLayer = apf::createNumbering(vMesh, "layer_num", apf::getLagrange(1) , 1); - apf::Numbering* nVID = apf::createNumbering(vMesh, "vid_num", apf::getLagrange(1) , 1); - - // process each layer - // layer = 1 corresponds to the surface nodes - // treat it out side of the loop - std::vector coords = readLayerCoordinates(1, t, s); - std::vector lastVs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, 1); - apf::number(nVID, v, 0, 0, i); - lastVs.push_back(v); - } - - PCU_ALWAYS_ASSERT(coords.size() == lastVs.size()); - - // now process the remaining layers - for (int layer = 2; layer < s+1; layer++) { - // read the coords first - std::vector coords = readLayerCoordinates(layer, t, s); - // create the new verts - std::vector vs; - for (std::size_t i = 0; i < coords.size(); i++) { - apf::MeshEntity* v = vMesh->createVertex(0, coords[i], apf::Vector3(0, 0, 0)); - apf::number(nLayer, v, 0, 0, layer); - apf::number(nVID, v, 0, 0, i); - vs.push_back(v); - } - - // create the prisms - for (std::size_t i = 0; i < faces.size(); i++) { - apf::MeshEntity* f = faces[i]; - apf::MeshEntity* dv[3]; - m->getDownward(f, 0, dv); - int id0 = apf::getNumber(n, dv[0], 0, 0); - int id1 = apf::getNumber(n, dv[1], 0, 0); - int id2 = apf::getNumber(n, dv[2], 0, 0); - - apf::MeshEntity* prismVs[6]; - - prismVs[0] = lastVs[id0]; - prismVs[1] = lastVs[id1]; - prismVs[2] = lastVs[id2]; - - prismVs[3] = vs[id0]; - prismVs[4] = vs[id1]; - prismVs[5] = vs[id2]; - - apf::buildElement(vMesh, 0, apf::Mesh::PRISM, prismVs); - } - // set the lastVs to Vs - for (std::size_t i = 0; i < vs.size(); i++) - lastVs[i] = vs[i]; - } - vMesh->acceptChanges(); - apf::deriveMdsModel(vMesh); - apf::verify(vMesh); - - // create the surface to strand map - it = vMesh->begin(0); - while ( (e = vMesh->iterate(it)) ) { - int layer = apf::getNumber(nLayer, e, 0, 0) - 1; - int vid = apf::getNumber(nVID, e, 0, 0); - PCU_ALWAYS_ASSERT((layer >= 0) && (layer < s)); - PCU_ALWAYS_ASSERT((vid >= 0) && (vid < (int)m->count(0))); - surfToStrandMap[vid][layer] = e; - } - vMesh->end(it); - - // return the volume mesh - return vMesh; -} - -void writeCre(CapstoneModule& cs, const std::string& fileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, fileName.c_str(), idmapping); -} - -void writeMdsMesh(apf::Mesh2* m, const char* name, const char* fieldName) -{ - // this would be the smb mesh with a field associated with it - apf::Mesh2* outMesh = apf::makeEmptyMdsMesh(gmi_load(".null"), 2, false); - - // add numbering to the surface mash - apf::Numbering* n = apf::createNumbering(m, "v_num", apf::getLagrange(1) , 1); - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(0); - int count = 0; - while ( (e = m->iterate(it)) ) { - apf::number(n, e, 0, 0, count); - count++; - } - m->end(it); - - - int numVerts = m->count(0); - apf::MeshEntity* vMap[numVerts]; - - apf::MeshEntity* v; - it = m->begin(0); - while( (v = m->iterate(it)) ) { - apf::Vector3 p; - m->getPoint(v, 0, p); - apf::MeshEntity* nv = outMesh->createVertex(0, p, apf::Vector3(0, 0, 0)); - int id = apf::getNumber(n, v, 0, 0); - vMap[id] = nv; - } - m->end(it); - - it = m->begin(2); - while( (e = m->iterate(it)) ) { - apf::MeshEntity* vs[3]; - m->getDownward(e, 0, vs); - apf::MeshEntity* triVs[3]; - triVs[0] = vMap[apf::getNumber(n, vs[0], 0, 0)]; - triVs[1] = vMap[apf::getNumber(n, vs[1], 0, 0)]; - triVs[2] = vMap[apf::getNumber(n, vs[2], 0, 0)]; - apf::buildElement(outMesh, 0, apf::Mesh::TRIANGLE, triVs); - } - m->end(it); - - outMesh->acceptChanges(); - apf::deriveMdsModel(outMesh); - outMesh->verify(); - apf::verify(outMesh); - - apf::destroyNumbering(n); - - // get the field on input mesh - apf::Field* inputField = m->findField(fieldName); - PCU_ALWAYS_ASSERT(inputField); - - // create the corresponding field for the out mesh - apf::Field* outputField = apf::createFieldOn(outMesh, "sizes", apf::SCALAR); - - it = m->begin(0); - apf::MeshIterator* outIt = outMesh->begin(0); - while( (v = m->iterate(it)) ) { - double value = apf::getScalar(inputField, v, 0); - apf::MeshEntity* outV = outMesh->iterate(outIt); - apf::setScalar(outputField, outV, 0, value); - } - m->end(it); - outMesh->end(outIt); - outMesh->writeNative(name); - - // clean up outMesh - outMesh->destroyNative(); - apf::destroyMesh(outMesh); -} diff --git a/capstone_clis/capStoneCurve.cc b/capstone_clis/capStoneCurve.cc deleted file mode 100644 index 8d525da10..000000000 --- a/capstone_clis/capStoneCurve.cc +++ /dev/null @@ -1,200 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// === includes for safe_mkdir === -#include -#include /*required for mode_t for mkdir on some systems*/ -#include /*using POSIX mkdir call for SMB "foo/" path*/ -#include /* for checking the error from mkdir */ - - -#include "CapstoneModule.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Analysis.h" -#include "CreateMG_Framework_Application.h" -#include "CreateMG_Framework_Attributes.h" -#include "CreateMG_Framework_Core.h" -#include "CreateMG_Framework_Geometry.h" -#include "CreateMG_Framework_Mesh.h" - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - -void safe_mkdir(const char* path) -{ - mode_t const mode = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH; - int err; - errno = 0; - err = mkdir(path, mode); - if (err != 0 && errno != EEXIST) - { - reel_fail("Err: could not create directory \"%s\"\n", path); - } -} - -void inspectFace(apf::Mesh2* m, apf::MeshEntity* e) -{ - apf::MeshEntity* down[12]; - int ndv = m->getDownward(e, 0, down); - int nde = m->getDownward(e, 1, down); - int ndf = m->getDownward(e, 2, down); - printf("downward v/e/f %d/%d/%d\n", ndv, nde, ndf); -} - -void inspectRegion(apf::Mesh2* m, apf::MeshEntity* e) -{ - apf::MeshEntity* down[12]; - int ndv = m->getDownward(e, 0, down); - int nde = m->getDownward(e, 1, down); - int ndf = m->getDownward(e, 2, down); - int ndr = m->getDownward(e, 3, down); - printf("downward v/e/f/r %d/%d/%d/%d\n", ndv, nde, ndf, ndr); -} - -void writeVtk(CapstoneModule& cs, const std::string& vtkFileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *vtkWriter = get_writer(ctx, "VTK File - Generic Writer"); - if (!vtkWriter) - error(HERE, ERR_GENERIC, "Could not find the VTK writer"); - - // This IO Tool re-indexes the vertices for the elements classified in the gregions set and writes out the VTK file - // the sets for edges and faces are empty so the mesh edge/faces classified over the egdes and faces wont get written explicitly - - /* std::vector vgregions; */ - /* vgregions.insert(vgregions.end(), gregions.begin(), gregions.end()); */ - /* vtkWriter->input().set_value("Regions", vgregions); */ - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - vtkWriter->write(ctx, gmodel, mmodels, vtkFileName.c_str(), idmapping); - - //cs.save_mesh_file(SystemTool::get_real_app_path().append("/booleaned.vtk"),ommodel); - -} - -int main(int argc, char** argv) -{ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - - lion_set_verbosity(1); - - gmi_cap_start(); - gmi_register_cap(); - - // load capstone mesh - - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("test", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(argv[1]); - - M_GModel gmodel = cs.load_files(filenames); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - - for(std::size_t i = 0; i < mmodels.size(); ++i) - { - M_MModel ammodel = mmodels[i]; - std::size_t numregs = 0; - MG_API_CALL(m, set_current_model(ammodel)); - MG_API_CALL(m, get_num_topos(TOPO_REGION, numregs)); - if(numregs > 0) - { - mmodel = ammodel; - break; - } - } - - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - /* writeVtk(cs, "before.vtk"); */ - apf::Mesh2* apfCapMesh = apf::createMesh(m, g); - - apf::Field* tf = apf::createFieldOn(apfCapMesh, "test_field", apf::VECTOR); - apf::MeshEntity* ent; - apf::MeshIterator* it = apfCapMesh->begin(0); - while( (ent = apfCapMesh->iterate(it)) ) { - apf::Vector3 p = ma::getPosition(apfCapMesh, ent); - double x = p[0]; - double y = p[1]; - double z = p[2]; - apf::Vector3 s(y, z, 2*x); - apf::setVector(tf, ent, 0, s); - } - apfCapMesh->end(it); - - apf::writeVtkFiles("before", apfCapMesh); - - int order = 2; - crv::BezierCurver bc(apfCapMesh, order, 0); - bc.run(); - - crv::writeCurvedVtuFiles(apfCapMesh, apf::Mesh::TRIANGLE, order + 2, "curved"); - crv::writeCurvedWireFrame(apfCapMesh, order + 10, "curved"); - - - M_MModel mmdl; - m->get_current_model(mmdl); - std::string info; - m->print_info(mmdl, info); - std::cout << info << std::endl; - - apf::writeVtkFiles("after", apfCapMesh); - /* writeVtk(cs, "after.vtk"); */ - - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} diff --git a/capstone_clis/capStoneIsoAdaptB737.cc b/capstone_clis/capStoneIsoAdaptB737.cc deleted file mode 100644 index e2cf07f00..000000000 --- a/capstone_clis/capStoneIsoAdaptB737.cc +++ /dev/null @@ -1,379 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "CapstoneModule.h" -/* #include "CreateMG_Framework_Core.h" */ -/* #include "CreateMG_Framework_Analysis.h" */ -/* #include "CreateMG_Framework_Application.h" */ -/* #include "CreateMG_Framework_Attributes.h" */ -/* #include "CreateMG_Framework_Core.h" */ -/* #include "CreateMG_Framework_Geometry.h" */ -/* #include "CreateMG_Framework_Mesh.h" */ - - -using namespace CreateMG; -using namespace CreateMG::Attribution; -using namespace CreateMG::Mesh; -using namespace CreateMG::Geometry; - -class B737 : public ma::IsotropicFunction -{ - public: - B737(ma::Mesh* m) - { - mesh = m; - } - virtual double getValue(ma::Entity* v) - { - // get current size - apf::Up up; - mesh->getUp(v, up); - double sum = 0; - double count = 0; - for (int i = 0; i < up.n; i++) { - ma::Entity* currentEdge = up.e[i]; - ma::Entity* vs[2]; - mesh->getDownward(currentEdge, 0, vs); - ma::Vector diff = ma::getPosition(mesh, vs[0]) - ma::getPosition(mesh, vs[1]); - sum += diff.getLength(); - count++; - } - double currentSize = sum/count; - - // compute the size multiplier - double y = 1.; - int tag = mesh->getModelTag(mesh->toModel(v)); - int type = mesh->getModelType(mesh->toModel(v)); - if (type == 2 && (tag == 66 || tag == 86)) - y = 3; - if (type == 1 && tag == 178) - y = 3; - // new size is multiplier * currentSize - return y*currentSize; - } - private: - ma::Mesh* mesh; -}; - - -void writeCre(CapstoneModule& cs, const std::string& vtkFileName) -{ - GeometryDatabaseInterface *gdbi = cs.get_geometry(); - MeshDatabaseInterface *mdbi = cs.get_mesh(); - AppContext *ctx = cs.get_context(); - - // Get the VTK writer. - Writer *creWriter = get_writer(ctx, "Create Native Writer"); - if (!creWriter) - error(HERE, ERR_GENERIC, "Could not find the CRE writer"); - - IdMapping idmapping; - std::vector mmodels; - M_GModel gmodel; - M_MModel mmodel; - gdbi->get_current_model(gmodel); - mdbi->get_current_model(mmodel); - mmodels.clear(); - mmodels.push_back(mmodel); - creWriter->write(ctx, gmodel, mmodels, vtkFileName.c_str(), idmapping); -} - -void gradeSize(apf::Mesh2* m, apf::Field* f, double beta) -{ - std::queue q; - apf::MeshEntity* e; - apf::MeshIterator* it = m->begin(1); - while( (e = m->iterate(it)) ) { - apf::MeshEntity* vs[2]; - m->getDownward(e, 0, vs); - double s0 = apf::getScalar(f, vs[0], 0); - double s1 = apf::getScalar(f, vs[1], 0); - double factor; - if (s0 > s1) factor = s0/s1; - else factor = s1/s0; - if (factor > beta) - q.push(e); - } - m->end(it); - - while(!q.empty()) { - apf::MeshEntity* e = q.front(); - q.pop(); - - apf::MeshEntity* vs[2]; - m->getDownward(e, 0, vs); - double s0 = apf::getScalar(f, vs[0], 0); - double s1 = apf::getScalar(f, vs[1], 0); - - apf::Up up; - if (s1 > beta*s0) { - s1 = beta*s0; - apf::setScalar(f, vs[1], 0, s1); - m->getUp(vs[1], up); - } - else if (s0 > beta*s1) { - s0 = beta*s1; - apf::setScalar(f, vs[0], 0, s0); - m->getUp(vs[0], up); - } - else - continue; - - for (int i = 0; i < up.n; i++) { - q.push(up.e[i]); - } - } -} - -static double triQuality(apf::Mesh2* m, apf::MeshEntity* tri) -{ - apf::Vector3 p0; - apf::Vector3 p1; - apf::Vector3 p2; - - apf::MeshEntity* vs[3]; - m->getDownward(tri, 0, vs); - - m->getPoint(vs[0], 0, p0); - m->getPoint(vs[1], 0, p1); - m->getPoint(vs[2], 0, p2); - - apf::Vector3 d1 = p1 - p0; - apf::Vector3 d2 = p2 - p0; - - apf::Vector3 c = apf::cross(d1, d2); - - double areaSquare = c*c / 4.; - - apf::MeshEntity* es[3]; - m->getDownward(tri, 1, es); - double sum = 0; - for (int i = 0; i < 3; i++) { - apf::MeshEntity* vs[2]; - m->getDownward(es[i], 0, vs); - apf::Vector3 p0; - apf::Vector3 p1; - m->getPoint(vs[0], 0, p0); - m->getPoint(vs[1], 0, p1); - sum += (p1-p0)*(p1-p0); - } - return 48*areaSquare/sum/sum; -} - - -void stats(apf::Mesh2* m, apf::Field* f, double desiredQ, double &minQ, int &nQ, int &nMinL, int &nMaxL) -{ - apf::MeshEntity* e; - apf::MeshIterator* it; - - nQ = nMinL = nMaxL = 0; - - minQ = 1000.0; - it = m->begin(2); - while ( (e = m->iterate(it)) ) { - double q = triQuality(m, e); - if (q < minQ) - minQ = q; - if (q < desiredQ) - nQ++; - } - m->end(it); - - it = m->begin(1); - while ( (e = m->iterate(it)) ) { - apf::MeshEntity* vs[2]; - m->getDownward(e, 0, vs); - apf::Vector3 p0; - apf::Vector3 p1; - m->getPoint(vs[0], 0, p0); - m->getPoint(vs[1], 0, p1); - double length = (p1-p0).getLength(); - double s0 = apf::getScalar(f, vs[0],0); - double s1 = apf::getScalar(f, vs[1],0); - double normalizedLength = 2*length/(s0+s1); - if (normalizedLength < 0.5) - nMinL++; - if (normalizedLength > 1.5) - nMaxL++; - } - m->end(it); -} - -int main(int argc, char** argv) -{ - /* INIT CALLS */ - MPI_Init(&argc, &argv); - PCU_Comm_Init(); - - if (argc < 2) { - if (PCU_Comm_Self() == 0) { - printf("USAGE: %s \n", argv[0]); - } - MPI_Finalize(); - exit(EXIT_FAILURE); - } - - const char* inFileName = argv[1]; - - gmi_cap_start(); - gmi_register_cap(); - - /* LOAD CAPSTONE MESH */ - // create an instance of the Capstone Module activating CREATE/CREATE/CREATE - // for the Geometry/Mesh/Attribution databases - const std::string gdbName("Geometry Database : SMLIB");// Switch Create with SMLIB for CAD - const std::string mdbName("Mesh Database : Create"); - const std::string adbName("Attribution Database : Create"); - - CapstoneModule cs("test", gdbName.c_str(), mdbName.c_str(), adbName.c_str()); - - GeometryDatabaseInterface *g = cs.get_geometry(); - MeshDatabaseInterface *m = cs.get_mesh(); - AppContext *c = cs.get_context(); - - PCU_ALWAYS_ASSERT(g); - PCU_ALWAYS_ASSERT(m); - PCU_ALWAYS_ASSERT(c); - - v_string filenames; - filenames.push_back(inFileName); - - M_GModel gmodel = cs.load_files(filenames); - - M_MModel mmodel; - // Pick the volume mesh model from associated mesh models to this geom model - std::vector mmodels; - MG_API_CALL(m, get_associated_mesh_models(gmodel, mmodels)); - PCU_ALWAYS_ASSERT(mmodels.size() == 1); - MG_API_CALL(m, set_current_model(mmodels[0])); - - /* SET THE ADJACENCIES */ - MG_API_CALL(m, set_adjacency_state(REGION2FACE| - REGION2EDGE| - REGION2VERTEX| - FACE2EDGE| - FACE2VERTEX)); - MG_API_CALL(m, set_reverse_states()); - MG_API_CALL(m, set_adjacency_scope(TOPO_EDGE, SCOPE_FULL)); - MG_API_CALL(m, set_adjacency_scope(TOPO_FACE, SCOPE_FULL)); - MG_API_CALL(m, compute_adjacency()); - - /* CONVERT THE MESH TO APF::MESH2 */ - apf::Mesh2* apfCapMesh = apf::createMesh(m, g); - - /* ADD A TEST FIELD TO THE MESH TO DEMONSTRATE SOLUTION TRANSFER */ - apf::Field* tf = apf::createFieldOn(apfCapMesh, "test_field", apf::VECTOR); - apf::MeshEntity* ent; - apf::MeshIterator* it = apfCapMesh->begin(0); - while( (ent = apfCapMesh->iterate(it)) ) { - apf::Vector3 p = ma::getPosition(apfCapMesh, ent); - double x = p[0]; - double y = p[1]; - double z = p[2]; - apf::Vector3 s(y, z, 2*x); - apf::setVector(tf, ent, 0, s); - } - apfCapMesh->end(it); - - /* WRITE THE BEFORE ADAPT MESH TO VTK USING APF VTK WRITER */ - apf::writeVtkFiles("before", apfCapMesh); - - - /* SIZE SETUP AND CALL ADAPT */ - - // setting the isotropic sizes - ma::IsotropicFunction* sf = new B737(apfCapMesh); - apf::Field* adaptSize = apf::createFieldOn(apfCapMesh, "adapt_size", apf::SCALAR); - apf::MeshEntity* v; - it = apfCapMesh->begin(0); - while( (v = apfCapMesh->iterate(it)) ) { - apf::setScalar(adaptSize, v, 0, sf->getValue(v)); - } - apfCapMesh->end(it); - apf::writeVtkFiles("size_before_grade", apfCapMesh); - - // grading the size field - gradeSize(apfCapMesh, adaptSize, 1.4); - apf::writeVtkFiles("size_after_grade", apfCapMesh); - - // adapt setup - ma::Input* in; - in = ma::makeAdvanced(ma::configure(apfCapMesh, adaptSize)); - in->shouldSnap = true; - in->shouldTransferParametric = true; - in->shouldFixShape = true; - in->shouldForceAdaptation = true; - if (apfCapMesh->getDimension() == 2) - in->goodQuality = 0.04; - else // 3D meshes - in->goodQuality = 0.027; // this is the mean-ratio cubed - in->maximumIterations = 10; - - ma::adaptVerbose(in, false); - - - /* CREATE A CONSTANT FIELD ON THE MESH WITH QUALITIES */ - apf::Field* qual = apf::createField(apfCapMesh, "quality", apf::SCALAR, apf::getConstant(2)); - apf::MeshEntity* f; - it = apfCapMesh->begin(2); - while ( (f = apfCapMesh->iterate(it)) ) { - double q = triQuality(apfCapMesh, f); - apf::setScalar(qual, f, 0, q); - } - apfCapMesh->end(it); - - - /* WRITE THE AFTER ADAPT MESH TO VTK USING APF VTK WRITER */ - apf::writeVtkFiles("after", apfCapMesh); - - - /* PRINTING OUT SOME STATES */ - size_t vCount = apfCapMesh->count(0); - size_t eCount = apfCapMesh->count(1); - size_t fCount = apfCapMesh->count(2); - - double minQ; - int nQ, nMinL, nMaxL; - - - stats(apfCapMesh, adaptSize, 0.2*0.2, minQ, nQ, nMinL, nMaxL); - - printf("vCount is %d\n", (int) vCount); - printf("eCount is %d\n", (int) eCount); - printf("fCount is %d\n", (int) fCount); - printf("minQ is %f\n", minQ); - printf("nQ is %d\n", nQ); - printf("nMinL is %d\n", nMinL); - printf("nMaxL is %d\n", nMaxL); - - /* PRINT ADAPTED MESH INFO */ - M_MModel mmdl; - m->get_current_model(mmdl); - std::string info; - m->print_info(mmdl, info); - - /* WRITE THE ADAPTED MESH IN NATIVE CREATE FORMAT */ - writeCre(cs, "after.cre"); - - /* CLEAN UPS */ - delete sf; - - /* EXIT CALLS */ - gmi_cap_stop(); - PCU_Comm_Free(); - MPI_Finalize(); -} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7c874b151..e0c303870 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -217,6 +217,9 @@ if(ENABLE_CAPSTONE) if(HAVE_CAPSTONE_SIZINGMETRICTOOL) util_exe_func(cap_smooth cap_smooth.cc) endif() + util_exe_func(cap2vtk cap2vtk.cc) + util_exe_func(capGeomTest capGeomTest.cc) + util_exe_func(capCheckParam capCheckParam.cc) endif() # send all the newly added utility executable targets diff --git a/capstone_clis/capStone2VTK.cc b/test/cap2vtk.cc similarity index 100% rename from capstone_clis/capStone2VTK.cc rename to test/cap2vtk.cc diff --git a/capstone_clis/capStoneCheckParametrization.cc b/test/capCheckParam.cc similarity index 100% rename from capstone_clis/capStoneCheckParametrization.cc rename to test/capCheckParam.cc diff --git a/capstone_clis/capStoneGeomTest.cc b/test/capGeomTest.cc similarity index 100% rename from capstone_clis/capStoneGeomTest.cc rename to test/capGeomTest.cc diff --git a/test/capVol.cc b/test/capVol.cc index a9149c585..8bd064686 100644 --- a/test/capVol.cc +++ b/test/capVol.cc @@ -24,7 +24,7 @@ using namespace CreateMG::Attribution; using namespace CreateMG::Mesh; using namespace CreateMG::Geometry; -#include "capStoneSizeFields.h" +#include "capVolSizeFields.h" namespace { diff --git a/capstone_clis/capStoneSizeFields.h b/test/capVolSizeFields.h similarity index 100% rename from capstone_clis/capStoneSizeFields.h rename to test/capVolSizeFields.h