From 062710808be29c2b2ae9ccd2ec26c4f430d3e839 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 17 Mar 2020 11:19:44 -0400 Subject: [PATCH 01/36] Adds required functionality for Bezier Fields --- crv/CMakeLists.txt | 1 + crv/crv.h | 3 +++ crv/crvBezierFields.cc | 22 ++++++++++++++++++++++ 3 files changed, 26 insertions(+) create mode 100644 crv/crvBezierFields.cc diff --git a/crv/CMakeLists.txt b/crv/CMakeLists.txt index 992589d23..d6ca17a57 100644 --- a/crv/CMakeLists.txt +++ b/crv/CMakeLists.txt @@ -20,6 +20,7 @@ set(SOURCES crvTables.cc crvQuality.cc crvVtk.cc + crvBezierFields.cc ) # Package headers diff --git a/crv/crv.h b/crv/crv.h index 60e8c2c1a..6232db17d 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -37,6 +37,9 @@ int getBlendingOrder(const int type); /** \brief count invalid elements of the mesh */ int countNumberInvalidElements(apf::Mesh2* m); +/** \brief converts filed f, which is interpolating to Bezier */ +void convertInterpolatingFieldToBezier(apf::Mesh2* m, apf::Field* f); + /** \brief Base Mesh curving object \details P is the order, S is the space dimension, different from the mesh dimension, used to distinguish between planar 2D diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc new file mode 100644 index 000000000..37527096f --- /dev/null +++ b/crv/crvBezierFields.cc @@ -0,0 +1,22 @@ +/* + * Copyright 2015 Scientific Computation Research Center + * + * This work is open source software, licensed under the terms of the + * BSD license as described in the LICENSE file in the top-level directory. + */ + +#include "crv.h" +#include "crvBezier.h" +#include "crvTables.h" +#include "crvQuality.h" +#include +#include + +namespace crv { + +void convertInterpolatingFieldToBezier(apf::Mesh2* m, apf::Field* f) +{ + // TODO: to be completed +} + +} From 61ae04bb4554d83e525f2a921109abe7c57af3ea Mon Sep 17 00:00:00 2001 From: Avinash Moharana Date: Wed, 18 Mar 2020 12:57:46 -0400 Subject: [PATCH 02/36] handle fields through Bezier shapes --- crv/CMakeLists.txt | 1 + crv/crv.h | 2 +- crv/crvBezierFields.cc | 103 ++++++++++++++++++++++++++++++++++++++++- 3 files changed, 104 insertions(+), 2 deletions(-) diff --git a/crv/CMakeLists.txt b/crv/CMakeLists.txt index d6ca17a57..235ee85f0 100644 --- a/crv/CMakeLists.txt +++ b/crv/CMakeLists.txt @@ -6,6 +6,7 @@ set(SOURCES crvBezier.cc crvBezierPoints.cc crvBezierShapes.cc + crvBezierFields.cc crvBlended.cc crvCurveMesh.cc crvElevation.cc diff --git a/crv/crv.h b/crv/crv.h index 6232db17d..6aa3e61e2 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -38,7 +38,7 @@ int getBlendingOrder(const int type); int countNumberInvalidElements(apf::Mesh2* m); /** \brief converts filed f, which is interpolating to Bezier */ -void convertInterpolatingFieldToBezier(apf::Mesh2* m, apf::Field* f); +void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f); /** \brief Base Mesh curving object \details P is the order, S is the space dimension, diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index 37527096f..d35391c6e 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -7,16 +7,117 @@ #include "crv.h" #include "crvBezier.h" +#include "crvShape.h" #include "crvTables.h" #include "crvQuality.h" #include #include +#include + namespace crv { -void convertInterpolatingFieldToBezier(apf::Mesh2* m, apf::Field* f) +static void convertVectorFieldInterpolationPoints(int n, int ne, + apf::NewArray& nodes, + apf::NewArray& c, + apf::NewArray& newNodes){ + + for(int i = 0; i < ne; ++i) + newNodes[i].zero(); + + for( int i = 0; i < ne; ++i) + for( int j = 0; j < n; ++j) + newNodes[i] += nodes[j]*c[i*n+j]; + +} + +static void convertScalarFieldInterpolationPoints(int n, int ne, + apf::NewArray& nodes, + apf::NewArray& c, + apf::NewArray& newNodes){ + + for(int i = 0; i < ne; ++i) + newNodes[i] = 0.; + + for( int i = 0; i < ne; ++i) + for( int j = 0; j < n; ++j) + newNodes[i] += nodes[j]*c[i*n+j]; + +} + +void convertInterpolationFieldPoints(apf::MeshEntity* e, + apf::Field* f, int n, int ne, apf::NewArray& c){ + + apf::NewArray l, b(ne); + apf::NewArray ls, bs(ne); + apf::Element* elem = + apf::createElement(f,e); + + if (apf::getValueType(f) == apf::VECTOR) { + apf::getVectorNodes(elem,l); + convertVectorFieldInterpolationPoints(n, ne, l, c, b); + for(int i = 0; i < ne; ++i) + apf::setVector(f, e, i, b[i]); + } + else if (apf::getValueType(f) == apf::SCALAR) { + apf::getScalarNodes(elem, ls); + convertScalarFieldInterpolationPoints(n, ne, ls, c, bs); + for(int i = 0; i < ne; ++i) + apf::setScalar(f, e, i, bs[i]); + } + else + printf("Field type not implemented\n"); + + apf::destroyElement(elem); +} + +void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) { // TODO: to be completed + apf::FieldShape * fs = apf::getShape(f); + int order = fs->getOrder(); + + int md = m_mesh->getDimension(); + int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]); + // go downward, and convert interpolating to control points + int startDim = md - (blendingOrder > 0); + + for(int d = startDim; d >= 1; --d){ + if(!fs->hasNodesIn(d)) continue; + int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); + int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); + apf::NewArray c; + getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[d],c); + apf::MeshEntity* e; + apf::MeshIterator* it = m_mesh->begin(d); + while ((e = m_mesh->iterate(it))){ + if(m_mesh->isOwned(e)) + convertInterpolationFieldPoints(e,f,n,ne,c); + } + m_mesh->end(it); + } + // if we have a full representation, we need to place internal nodes on + // triangles and tetrahedra + for(int d = 2; d <= md; ++d){ + if(!fs->hasNodesIn(d) || + getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; + int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); + int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); + apf::NewArray c; + getInternalBezierTransformationCoefficients(m_mesh,order,1, + apf::Mesh::simplexTypes[d],c); + apf::MeshEntity* e; + apf::MeshIterator* it = m_mesh->begin(d); + while ((e = m_mesh->iterate(it))){ + if(!isBoundaryEntity(m_mesh,e) && m_mesh->isOwned(e)) + convertInterpolationFieldPoints(e,f,n-ne,ne,c); + } + m_mesh->end(it); + } + + //synchronize(); + } } From 1d6ab96334d6321e74870b7b7926b0bcc3ee03d9 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Mon, 23 Mar 2020 10:58:42 -0400 Subject: [PATCH 03/36] Takes care of blending when doing Bezier Fields Note that blending is used internally in Bezier Shapes. More specifically when a BezierCurver object is set with a non-zero third argument, Bezier shapes for tets and triangles might behave differently. Therefore, we needed to make sure the blending orders are set to zero at the beginning of the convertInterpolationgFieldsToBezier routine and then reverted back to the original values at the end. There still seems to be necessary to construct a new BezierCurver object with last argument 0 for the bezier fields to work properly. --- crv/crvBezierFields.cc | 39 ++++++++++++++++++--------------------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index d35391c6e..c2a3d93f3 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -10,6 +10,7 @@ #include "crvShape.h" #include "crvTables.h" #include "crvQuality.h" +#include #include #include @@ -50,8 +51,10 @@ void convertInterpolationFieldPoints(apf::MeshEntity* e, apf::NewArray l, b(ne); apf::NewArray ls, bs(ne); + apf::MeshElement* me = + apf::createMeshElement(f->getMesh(),e); apf::Element* elem = - apf::createElement(f,e); + apf::createElement(f,me); if (apf::getValueType(f) == apf::VECTOR) { apf::getVectorNodes(elem,l); @@ -67,20 +70,29 @@ void convertInterpolationFieldPoints(apf::MeshEntity* e, } else printf("Field type not implemented\n"); - apf::destroyElement(elem); + apf::destroyMeshElement(me); } void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) { // TODO: to be completed + // Zero out the blending orders (since they are used internally + // in Bezier shapes (triangles and tets) but save the old ones to + // revert back to at the end of this (so the rest of Bezier related + // procedures work as intended) + int oldBlendings[apf::Mesh::TYPES]; + for (int i = 0; i < apf::Mesh::TYPES; i++) + oldBlendings[i] = getBlendingOrder(i); + setBlendingOrder(apf::Mesh::TYPES, 0); + apf::FieldShape * fs = apf::getShape(f); int order = fs->getOrder(); int md = m_mesh->getDimension(); int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]); // go downward, and convert interpolating to control points - int startDim = md - (blendingOrder > 0); + int startDim = md - (blendingOrder > 0); for(int d = startDim; d >= 1; --d){ if(!fs->hasNodesIn(d)) continue; @@ -97,24 +109,9 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) } m_mesh->end(it); } - // if we have a full representation, we need to place internal nodes on - // triangles and tetrahedra - for(int d = 2; d <= md; ++d){ - if(!fs->hasNodesIn(d) || - getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; - int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); - int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); - apf::NewArray c; - getInternalBezierTransformationCoefficients(m_mesh,order,1, - apf::Mesh::simplexTypes[d],c); - apf::MeshEntity* e; - apf::MeshIterator* it = m_mesh->begin(d); - while ((e = m_mesh->iterate(it))){ - if(!isBoundaryEntity(m_mesh,e) && m_mesh->isOwned(e)) - convertInterpolationFieldPoints(e,f,n-ne,ne,c); - } - m_mesh->end(it); - } + + for (int i = 0; i < apf::Mesh::TYPES; i++) + setBlendingOrder(i, oldBlendings[i]); //synchronize(); From 775b74368739f9c080754e26e2bf51fbcb79e025 Mon Sep 17 00:00:00 2001 From: Avinash Moharana Date: Tue, 24 Mar 2020 12:19:24 -0400 Subject: [PATCH 04/36] bezier representation of solution fields --- crv/crvBezierFields.cc | 40 +++++++++++++++------------------------- crv/crvCurveMesh.cc | 1 + 2 files changed, 16 insertions(+), 25 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index d35391c6e..1aca7e1df 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -7,11 +7,14 @@ #include "crv.h" #include "crvBezier.h" +#include "crvBezierShapes.h" +#include "crvMath.h" #include "crvShape.h" #include "crvTables.h" #include "crvQuality.h" #include #include +#include #include @@ -28,7 +31,6 @@ static void convertVectorFieldInterpolationPoints(int n, int ne, for( int i = 0; i < ne; ++i) for( int j = 0; j < n; ++j) newNodes[i] += nodes[j]*c[i*n+j]; - } static void convertScalarFieldInterpolationPoints(int n, int ne, @@ -42,11 +44,11 @@ static void convertScalarFieldInterpolationPoints(int n, int ne, for( int i = 0; i < ne; ++i) for( int j = 0; j < n; ++j) newNodes[i] += nodes[j]*c[i*n+j]; - } void convertInterpolationFieldPoints(apf::MeshEntity* e, - apf::Field* f, int n, int ne, apf::NewArray& c){ + apf::Field* f, + int n, int ne, apf::NewArray& c){ apf::NewArray l, b(ne); apf::NewArray ls, bs(ne); @@ -77,12 +79,18 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) apf::FieldShape * fs = apf::getShape(f); int order = fs->getOrder(); + //apf::Field* fnew = createField( + // m_mesh, "copy_field", apf::getValueType(f), crv::getBezier(order)); + //transferFields(m_mesh, f, fnew); + int md = m_mesh->getDimension(); - int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]); + //int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]); + //blendingOrder = 0; // go downward, and convert interpolating to control points - int startDim = md - (blendingOrder > 0); + //int startDim = md - (blendingOrder > 0); - for(int d = startDim; d >= 1; --d){ + //for(int d = startDim; d >= 1; --d){ + for(int d = md; d >= 1; d--){ if(!fs->hasNodesIn(d)) continue; int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); @@ -97,26 +105,8 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) } m_mesh->end(it); } - // if we have a full representation, we need to place internal nodes on - // triangles and tetrahedra - for(int d = 2; d <= md; ++d){ - if(!fs->hasNodesIn(d) || - getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; - int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); - int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); - apf::NewArray c; - getInternalBezierTransformationCoefficients(m_mesh,order,1, - apf::Mesh::simplexTypes[d],c); - apf::MeshEntity* e; - apf::MeshIterator* it = m_mesh->begin(d); - while ((e = m_mesh->iterate(it))){ - if(!isBoundaryEntity(m_mesh,e) && m_mesh->isOwned(e)) - convertInterpolationFieldPoints(e,f,n-ne,ne,c); - } - m_mesh->end(it); - } - //synchronize(); + apf::synchronize(f); } diff --git a/crv/crvCurveMesh.cc b/crv/crvCurveMesh.cc index 5cc39bbfb..966706b30 100644 --- a/crv/crvCurveMesh.cc +++ b/crv/crvCurveMesh.cc @@ -150,6 +150,7 @@ void BezierCurver::convertInterpolatingToBezier() } // if we have a full representation, we need to place internal nodes on // triangles and tetrahedra + for(int d = 2; d <= md; ++d){ if(!fs->hasNodesIn(d) || getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; From 2982d4ef958e755b3fcc75039cbc93ffe96b43d3 Mon Sep 17 00:00:00 2001 From: Avinash Date: Tue, 7 Apr 2020 15:43:34 -0400 Subject: [PATCH 05/36] testing bezier field solution transfer. --- crv/crv.h | 6 +- crv/crvBezierFieldTransfer.cc | 117 ++++++++++++++++++++++++++++++++++ crv/crvBezierFieldTransfer.h | 43 +++++++++++++ 3 files changed, 165 insertions(+), 1 deletion(-) create mode 100644 crv/crvBezierFieldTransfer.cc create mode 100644 crv/crvBezierFieldTransfer.h diff --git a/crv/crv.h b/crv/crv.h index 6aa3e61e2..00f166912 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -37,7 +37,11 @@ int getBlendingOrder(const int type); /** \brief count invalid elements of the mesh */ int countNumberInvalidElements(apf::Mesh2* m); -/** \brief converts filed f, which is interpolating to Bezier */ +/** \ brief converts field f to Bezier entity wise */ +void convertInterpolationFieldPoints(apf::MeshEntity* e, + apf::Field* f, int n, int ne, apf::NewArray &c); + +/** \brief converts field f, which is interpolating to Bezier */ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f); /** \brief Base Mesh curving object diff --git a/crv/crvBezierFieldTransfer.cc b/crv/crvBezierFieldTransfer.cc new file mode 100644 index 000000000..3bb8ee3e4 --- /dev/null +++ b/crv/crvBezierFieldTransfer.cc @@ -0,0 +1,117 @@ +/* +1;3409;0c * Copyright 2015 Scientific Computation Research Center + * + * This work is open source software, licensed under the terms of the + * BSD license as described in the LICENSE file in the top-level directory. + */ +#include + +#include "crvAdapt.h" +#include "crvBezier.h" +#include "crvBezierShapes.h" +#include "crvMath.h" +#include "crvQuality.h" +#include "crvShape.h" +#include "crvTables.h" +#include +#include +#include +#include +#include +#include + +namespace crv { +class LinearTransfer; +class CavityTransfer; + +CrvBezierFieldTransfer::~CrvBezierFieldTransfer() +{ +} + +void CrvBezierFieldTransfer::onVertex( + apf::MeshElement*, + Vector const&, + Entity*) +{ +} + +void CrvBezierFieldTransfer::onRefine( + Entity*, + EntityArray&) +{ +} + +void CrvBezierFieldTransfer::onCavity( + EntityArray&, + EntityArray&) +{ +} + +class HigherOrderTransfer : public CrvBezierField +{ + public: + LinearTransfer verts; + CavityTransfer others; + HighOrderTransfer(apf::Field* f): + verts(f),others(f) + { + } + virtual bool hasNodesOn(int dimension) + { + return others.hasNodesOn(dimension); + } + virtual void onVertex( + apf::MeshElement* parent, + Vector const& xi, + Entity* vert) + { + verts.onVertex(parent,xi,vert); + } + virtual void onRefine( + Entity* parent, + EntityArray& newEntities) + { + others.onRefine(parent,newEntities); + } + virtual void onCavity( + EntityArray& oldElements, + EntityArray& newEntities) + { + others.onCavity(oldElements,newEntities); + } + virtual void convertToBezierFields( + int dimension, + EntityArray& newEntities) + { + //apf::FieldShape *fs = others.shape; + //std::string name = fs->getName(); + //if (name == std::string("Bezier")) { + int order = shape->getOrder(); + int n = shape->getEntityShape( + apf::Mesh::simplexTypes[dimension])->countNodes(); + int ne = shape->countNodesOn( + apf::Mesh::simplexTypes[dimension]); + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[dimension], c); + + for( size_t i = 0; i < newEntities.getSize(); i++) { + crv::convertInterpolationFieldPoints(newEntities[i], + field, n, ne, c); + } + } +}; + +CrvBezierFieldTransfer::CrvBezierFieldTransfer(Mesh* m) +{ + for (int i = 0; i < m->countFields(); ++i) + { + apf::Field* f = m->getField(i); + printf("Fields considered, %s\n", apf::getName(f)); + this->add(createFieldTransfer(f)); + } +} + + +} + diff --git a/crv/crvBezierFieldTransfer.h b/crv/crvBezierFieldTransfer.h new file mode 100644 index 000000000..d0f381456 --- /dev/null +++ b/crv/crvBezierFieldTransfer.h @@ -0,0 +1,43 @@ +/* + * Copyright 2015 Scientific Computation Research Center + * + * This work is open source software, licensed under the terms of the + * BSD license as described in the LICENSE file in the top-level directory. + */ + +#ifndef CRVBEZIERSHAPES_H +#define CRVBEZIERSHAPES_H + +/** \file crvBezierShapes.h + * \brief main file for bezier shape functions */ + +namespace crv { +class CrvBezierFieldTransfer +{ + public: + virtual ~CrvBezierFieldTransfer(); + + virtual bool hasNodesOn(int dimension) = 0; + + virtual void onVertex( + apf::MeshElement* parent, + Vector const& xi, + Entity* vert); + + virtual void onRefine( + Entity* parent, + EntityArray& newEntities); + + virtual void onCavity( + EntityArray& oldElements, + EntityArray& newEntities); + + virtual void convertToBezierFields( + int dimension, + EntityArray &newEntities); +} + +CrvBezierFieldTransfer* createFieldTransfer(f); + + +#endif From 47fb3b6962df9aae851a78df623048ce56297419 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 7 Apr 2020 15:51:41 -0400 Subject: [PATCH 06/36] Renames files FieldTransfer to SolutionTransfer --- crv/{crvBezierFieldTransfer.cc => crvBezierSolutionTransfer.cc} | 0 crv/{crvBezierFieldTransfer.h => crvBezierSolutionTransfer.h} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename crv/{crvBezierFieldTransfer.cc => crvBezierSolutionTransfer.cc} (100%) rename crv/{crvBezierFieldTransfer.h => crvBezierSolutionTransfer.h} (100%) diff --git a/crv/crvBezierFieldTransfer.cc b/crv/crvBezierSolutionTransfer.cc similarity index 100% rename from crv/crvBezierFieldTransfer.cc rename to crv/crvBezierSolutionTransfer.cc diff --git a/crv/crvBezierFieldTransfer.h b/crv/crvBezierSolutionTransfer.h similarity index 100% rename from crv/crvBezierFieldTransfer.h rename to crv/crvBezierSolutionTransfer.h From 2228cfcb0f4957b378e5fa33ba4ebae625ff03ac Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Tue, 7 Apr 2020 16:36:16 -0400 Subject: [PATCH 07/36] Adds setBezierSolutionTransfers for clients --- crv/crv.h | 5 +++ crv/crvBezierSolutionTransfer.cc | 59 ++++++++++++-------------------- crv/crvBezierSolutionTransfer.h | 10 +++--- 3 files changed, 31 insertions(+), 43 deletions(-) diff --git a/crv/crv.h b/crv/crv.h index 00f166912..e8f0d1fe5 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -11,6 +11,7 @@ #include "apfMesh2.h" #include "apfShape.h" #include +#include #include #include #include @@ -210,6 +211,10 @@ void writeInterpolationPointVtuFiles(apf::Mesh* m, const char* prefix); int getTriNodeIndex(int P, int i, int j); int getTetNodeIndex(int P, int i, int j, int k); +/** \brief adds bezier solution transfers */ +ma::SolutionTransfer* setBezierSolutionTransfers( + const std::vector& fields); + /** \brief crv fail function */ void fail(const char* why) __attribute__((noreturn)); diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 3bb8ee3e4..f265f338a 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -13,6 +13,7 @@ #include "crvQuality.h" #include "crvShape.h" #include "crvTables.h" +#include "crvBezierSolutionTransfer.h" #include #include #include @@ -21,38 +22,15 @@ #include namespace crv { -class LinearTransfer; -class CavityTransfer; +class ma::LinearTransfer; +class ma::CavityTransfer; -CrvBezierFieldTransfer::~CrvBezierFieldTransfer() -{ -} - -void CrvBezierFieldTransfer::onVertex( - apf::MeshElement*, - Vector const&, - Entity*) -{ -} - -void CrvBezierFieldTransfer::onRefine( - Entity*, - EntityArray&) -{ -} - -void CrvBezierFieldTransfer::onCavity( - EntityArray&, - EntityArray&) -{ -} - -class HigherOrderTransfer : public CrvBezierField +class CrvBezierSolutionTransfer : public ma::SolutionTransfer { public: - LinearTransfer verts; - CavityTransfer others; - HighOrderTransfer(apf::Field* f): + ma::LinearTransfer verts; + ma::CavityTransfer others; + CrvBezierSolutionTransfer(apf::Field* f): verts(f),others(f) { } @@ -79,13 +57,15 @@ class HigherOrderTransfer : public CrvBezierField { others.onCavity(oldElements,newEntities); } - virtual void convertToBezierFields( + protected: + void convertToBezierFields( int dimension, EntityArray& newEntities) { //apf::FieldShape *fs = others.shape; //std::string name = fs->getName(); //if (name == std::string("Bezier")) { + int td = apf::Mesh::typeDimension[] int order = shape->getOrder(); int n = shape->getEntityShape( apf::Mesh::simplexTypes[dimension])->countNodes(); @@ -102,16 +82,19 @@ class HigherOrderTransfer : public CrvBezierField } }; -CrvBezierFieldTransfer::CrvBezierFieldTransfer(Mesh* m) +static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f) { - for (int i = 0; i < m->countFields(); ++i) - { - apf::Field* f = m->getField(i); - printf("Fields considered, %s\n", apf::getName(f)); - this->add(createFieldTransfer(f)); - } + return new CrvBezierSolutionTransfer(f); } - +ma::SolutionTransfer* setBezierSolutionTransfers( + const std::vector& fields) +{ + ma::SolutionTransfer* st = new ma::SolutionTransfers(); + for (std::size_t i = 0; i < fields.size(); i++) { + st->add(createBezierSolutionTransfer(fields[i])); + } + return st; } +} diff --git a/crv/crvBezierSolutionTransfer.h b/crv/crvBezierSolutionTransfer.h index d0f381456..9bfd11d75 100644 --- a/crv/crvBezierSolutionTransfer.h +++ b/crv/crvBezierSolutionTransfer.h @@ -5,17 +5,17 @@ * BSD license as described in the LICENSE file in the top-level directory. */ -#ifndef CRVBEZIERSHAPES_H -#define CRVBEZIERSHAPES_H +#ifndef CRVBEZIERSOLUTIONTRANSFER_H +#define CRVBEZIERSOLUTIONTRANSFER_H /** \file crvBezierShapes.h * \brief main file for bezier shape functions */ namespace crv { -class CrvBezierFieldTransfer +class CrvBezierSolutionTransfer { public: - virtual ~CrvBezierFieldTransfer(); + virtual ~CrvBezierSolutionTransfer(); virtual bool hasNodesOn(int dimension) = 0; @@ -37,7 +37,7 @@ class CrvBezierFieldTransfer EntityArray &newEntities); } -CrvBezierFieldTransfer* createFieldTransfer(f); +CrvBezierSolutionTransfer* createFieldTransfer(f); #endif From f71dfeee69d8b2758af5f0e48758074c12def4da Mon Sep 17 00:00:00 2001 From: Avinash Date: Wed, 8 Apr 2020 13:15:59 -0400 Subject: [PATCH 08/36] bezier field representation after Refine. --- crv/crvBezierFields.cc | 6 +++--- crv/crvBezierSolutionTransfer.cc | 21 ++++++++++++++++++++- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index d35054279..06cabcfd6 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -83,9 +83,9 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) // in Bezier shapes (triangles and tets) but save the old ones to // revert back to at the end of this (so the rest of Bezier related // procedures work as intended) - int oldBlendings[apf::Mesh::TYPES]; - for (int i = 0; i < apf::Mesh::TYPES; i++) - oldBlendings[i] = getBlendingOrder(i); + //int oldBlendings[apf::Mesh::TYPES]; + //for (int i = 0; i < apf::Mesh::TYPES; i++) + // oldBlendings[i] = getBlendingOrder(i); setBlendingOrder(apf::Mesh::TYPES, 0); apf::FieldShape * fs = apf::getShape(f); diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index f265f338a..5681c18db 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -6,6 +6,7 @@ */ #include +#include "crv.h" #include "crvAdapt.h" #include "crvBezier.h" #include "crvBezierShapes.h" @@ -13,7 +14,6 @@ #include "crvQuality.h" #include "crvShape.h" #include "crvTables.h" -#include "crvBezierSolutionTransfer.h" #include #include #include @@ -50,6 +50,23 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer EntityArray& newEntities) { others.onRefine(parent,newEntities); + apf::FieldShape *shape = others.shape; + apf::Mesh *m = others.mesh; + int type = m->getType(newEntitites[0]); + int td = apf::Mesh::typeDimension[type]; + int order = shape->getOrder(); + int n = shape->getEntityShape( + apf::Mesh::simplexTypes[td])->countNodes(); + int ne = shape->countNodesOn( + apf::Mesh::simplexTypes[td]); + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[td], c); + + for( size_t i = 0; i < newEntities.getSize(); i++) { + crv::convertInterpolationFieldPoints(newEntities[i], + field, n, ne, c); + } } virtual void onCavity( EntityArray& oldElements, @@ -57,6 +74,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer { others.onCavity(oldElements,newEntities); } + /* protected: void convertToBezierFields( int dimension, @@ -80,6 +98,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer field, n, ne, c); } } + */ }; static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f) From 9a09525f0631884290ddfb1102d2827dc0390ae8 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Wed, 8 Apr 2020 14:54:26 -0400 Subject: [PATCH 09/36] Separates implementation of maSolutionTransfer --- ma/maSolutionTransfer.cc | 235 +++++++++++++++++---------------------- ma/maSolutionTransfer.h | 61 ++++++++++ 2 files changed, 166 insertions(+), 130 deletions(-) diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index 89b564e82..e7ae9f403 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -10,9 +10,6 @@ #include "maSolutionTransfer.h" #include "maAffine.h" #include "maMap.h" -#include -#include -#include namespace ma { @@ -63,142 +60,120 @@ int SolutionTransfer::getTransferDimension() return transferDimension; } -class FieldTransfer : public SolutionTransfer + +FieldTransfer::FieldTransfer(apf::Field* f) { - public: - FieldTransfer(apf::Field* f) - { - field = f; - mesh = apf::getMesh(f); - shape = apf::getShape(f); - value.allocate(apf::countComponents(f)); - } - /* hmm... in vs. on ... probably the ma:: signature - should change, it has the least users */ - virtual bool hasNodesOn(int dimension) - { - return shape->hasNodesIn(dimension); - } - apf::Field* field; - apf::Mesh* mesh; - apf::FieldShape* shape; - apf::NewArray value; -}; + field = f; + mesh = apf::getMesh(f); + shape = apf::getShape(f); + value.allocate(apf::countComponents(f)); +} -class LinearTransfer : public FieldTransfer +bool FieldTransfer::hasNodesOn(int dimension) { - public: - LinearTransfer(apf::Field* f): - FieldTransfer(f) - { - } - virtual void onVertex( - apf::MeshElement* parent, - Vector const& xi, - Entity* vert) - { - apf::Element* e = apf::createElement(field,parent); - apf::getComponents(e,xi,&(value[0])); - apf::setComponents(field,vert,0,&(value[0])); - apf::destroyElement(e); - } -}; + return shape->hasNodesIn(dimension); +} -class CavityTransfer : public FieldTransfer +void LinearTransfer::onVertex( + apf::MeshElement* parent, + Vector const& xi, + Entity* vert) { - public: - CavityTransfer(apf::Field* f): - FieldTransfer(f) - { - minDim = getMinimumDimension(getShape(f)); - } - void transferToNodeIn( - apf::Element* elem, - apf::Node const& node, - Vector const& elemXi) - { - apf::getComponents(elem,elemXi,&(value[0])); - apf::setComponents(field,node.entity,node.node,&(value[0])); - } - int getBestElement( - int n, - apf::Element** elems, - Affine* elemInvMaps, - Vector const& point, - Vector& bestXi) - { - double bestValue = -DBL_MAX; - int bestI = 0; - for (int i = 0; i < n; ++i) - { - Vector xi = elemInvMaps[i] * point; - double value = getInsideness(mesh,apf::getMeshEntity(elems[i]),xi); - if (value > bestValue) - { - bestValue = value; - bestI = i; - bestXi = xi; - } - } - return bestI; - } - void transferToNode( - int n, - apf::Element** elems, - Affine* elemInvMaps, - apf::Node const& node) - { - Vector xi; - shape->getNodeXi(mesh->getType(node.entity),node.node,xi); - Affine childMap = getMap(mesh,node.entity); - Vector point = childMap * xi; - Vector elemXi; - int i = getBestElement(n,elems,elemInvMaps,point,elemXi); - transferToNodeIn(elems[i],node,elemXi); - } - void transfer( - int n, - Entity** cavity, - EntityArray& newEntities) - { - if (getDimension(mesh, cavity[0]) < minDim) - return; - apf::NewArray elems(n); - for (int i = 0; i < n; ++i) - elems[i] = apf::createElement(field,cavity[i]); - apf::NewArray elemInvMaps(n); - for (int i = 0; i < n; ++i) - elemInvMaps[i] = invert(getMap(mesh,cavity[i])); - for (size_t i = 0; i < newEntities.getSize(); ++i) - { - int type = mesh->getType(newEntities[i]); - if (type == apf::Mesh::VERTEX) - continue; //vertices will have been handled specially beforehand - int nnodes = shape->countNodesOn(type); - for (int j = 0; j < nnodes; ++j) - { - apf::Node node(newEntities[i],j); - transferToNode(n,&(elems[0]),&(elemInvMaps[0]),node); - } - } - for (int i = 0; i < n; ++i) - apf::destroyElement(elems[i]); - } - virtual void onRefine( - Entity* parent, - EntityArray& newEntities) + apf::Element* e = apf::createElement(field,parent); + apf::getComponents(e,xi,&(value[0])); + apf::setComponents(field,vert,0,&(value[0])); + apf::destroyElement(e); +} + +CavityTransfer::CavityTransfer(apf::Field* f): + FieldTransfer(f) +{ + minDim = getMinimumDimension(getShape(f)); +} +void CavityTransfer::transferToNodeIn( + apf::Element* elem, + apf::Node const& node, + Vector const& elemXi) +{ + apf::getComponents(elem,elemXi,&(value[0])); + apf::setComponents(field,node.entity,node.node,&(value[0])); +} +int CavityTransfer::getBestElement( + int n, + apf::Element** elems, + Affine* elemInvMaps, + Vector const& point, + Vector& bestXi) +{ + double bestValue = -DBL_MAX; + int bestI = 0; + for (int i = 0; i < n; ++i) + { + Vector xi = elemInvMaps[i] * point; + double value = getInsideness(mesh,apf::getMeshEntity(elems[i]),xi); + if (value > bestValue) { - transfer(1,&parent,newEntities); + bestValue = value; + bestI = i; + bestXi = xi; } - virtual void onCavity( - EntityArray& oldElements, - EntityArray& newEntities) + } + return bestI; +} +void CavityTransfer::transferToNode( + int n, + apf::Element** elems, + Affine* elemInvMaps, + apf::Node const& node) +{ + Vector xi; + shape->getNodeXi(mesh->getType(node.entity),node.node,xi); + Affine childMap = getMap(mesh,node.entity); + Vector point = childMap * xi; + Vector elemXi; + int i = getBestElement(n,elems,elemInvMaps,point,elemXi); + transferToNodeIn(elems[i],node,elemXi); +} +void CavityTransfer::transfer( + int n, + Entity** cavity, + EntityArray& newEntities) +{ + if (getDimension(mesh, cavity[0]) < minDim) + return; + apf::NewArray elems(n); + for (int i = 0; i < n; ++i) + elems[i] = apf::createElement(field,cavity[i]); + apf::NewArray elemInvMaps(n); + for (int i = 0; i < n; ++i) + elemInvMaps[i] = invert(getMap(mesh,cavity[i])); + for (size_t i = 0; i < newEntities.getSize(); ++i) + { + int type = mesh->getType(newEntities[i]); + if (type == apf::Mesh::VERTEX) + continue; //vertices will have been handled specially beforehand + int nnodes = shape->countNodesOn(type); + for (int j = 0; j < nnodes; ++j) { - transfer(oldElements.getSize(),&(oldElements[0]),newEntities); + apf::Node node(newEntities[i],j); + transferToNode(n,&(elems[0]),&(elemInvMaps[0]),node); } - private: - int minDim; -}; + } + for (int i = 0; i < n; ++i) + apf::destroyElement(elems[i]); +} +void CavityTransfer::onRefine( + Entity* parent, + EntityArray& newEntities) +{ + transfer(1,&parent,newEntities); +} +void CavityTransfer::onCavity( + EntityArray& oldElements, + EntityArray& newEntities) +{ + transfer(oldElements.getSize(),&(oldElements[0]),newEntities); +} /* hmm... could use multiple inheritance here, but that creates a "diamond problem": diff --git a/ma/maSolutionTransfer.h b/ma/maSolutionTransfer.h index a80fe90e5..1d59e65b4 100644 --- a/ma/maSolutionTransfer.h +++ b/ma/maSolutionTransfer.h @@ -14,10 +14,15 @@ \brief MeshAdapt solution transfer interface */ #include +#include +#include +#include #include "maMesh.h" namespace ma { +class Affine; + /** \brief user-defined solution transfer base \details one of these classes should be defined for each field the user wants transferred during @@ -73,6 +78,62 @@ class SolutionTransfer int getTransferDimension(); }; +class FieldTransfer : public SolutionTransfer +{ + public: + FieldTransfer(apf::Field* f); + virtual bool hasNodesOn(int dimension); + apf::Field* field; + apf::Mesh* mesh; + apf::FieldShape* shape; + apf::NewArray value; +}; + +class LinearTransfer : public FieldTransfer +{ + public: + LinearTransfer(apf::Field* f): + FieldTransfer(f) {} + virtual void onVertex( + apf::MeshElement* parent, + Vector const& xi, + Entity* vert); +}; + +class CavityTransfer : public FieldTransfer +{ + public: + CavityTransfer(apf::Field* f); + void transferToNodeIn( + apf::Element* elem, + apf::Node const& node, + Vector const& elemXi); + int getBestElement( + int n, + apf::Element** elems, + Affine* elemInvMaps, + Vector const& point, + Vector& bestXi); + void transferToNode( + int n, + apf::Element** elems, + Affine* elemInvMaps, + apf::Node const& node); + void transfer( + int n, + Entity** cavity, + EntityArray& newEntities); + virtual void onRefine( + Entity* parent, + EntityArray& newEntities); + virtual void onCavity( + EntityArray& oldElements, + EntityArray& newEntities); + private: + int minDim; +}; + + /** \brief Creates a default solution transfer object for a field \details MeshAdapt has good algorithms for transferring nodal fields as well as using the Voronoi system for transferring From ac87f00f230010ccbccd68e101198b4354f3e208 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Wed, 8 Apr 2020 14:55:22 -0400 Subject: [PATCH 10/36] Adds crvBezierSolutionTransfer to CMake --- crv/CMakeLists.txt | 1 + crv/crvBezierSolutionTransfer.cc | 21 ++++++++++----------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/crv/CMakeLists.txt b/crv/CMakeLists.txt index 235ee85f0..86f650ae0 100644 --- a/crv/CMakeLists.txt +++ b/crv/CMakeLists.txt @@ -7,6 +7,7 @@ set(SOURCES crvBezierPoints.cc crvBezierShapes.cc crvBezierFields.cc + crvBezierSolutionTransfer.cc crvBlended.cc crvCurveMesh.cc crvElevation.cc diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 5681c18db..a141c1878 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -14,6 +14,7 @@ #include "crvQuality.h" #include "crvShape.h" #include "crvTables.h" +#include #include #include #include @@ -22,8 +23,6 @@ #include namespace crv { -class ma::LinearTransfer; -class ma::CavityTransfer; class CrvBezierSolutionTransfer : public ma::SolutionTransfer { @@ -40,19 +39,19 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer } virtual void onVertex( apf::MeshElement* parent, - Vector const& xi, - Entity* vert) + ma::Vector const& xi, + ma::Entity* vert) { verts.onVertex(parent,xi,vert); } virtual void onRefine( - Entity* parent, - EntityArray& newEntities) + ma::Entity* parent, + ma::EntityArray& newEntities) { others.onRefine(parent,newEntities); apf::FieldShape *shape = others.shape; apf::Mesh *m = others.mesh; - int type = m->getType(newEntitites[0]); + int type = m->getType(newEntities[0]); int td = apf::Mesh::typeDimension[type]; int order = shape->getOrder(); int n = shape->getEntityShape( @@ -65,12 +64,12 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer for( size_t i = 0; i < newEntities.getSize(); i++) { crv::convertInterpolationFieldPoints(newEntities[i], - field, n, ne, c); + others.field, n, ne, c); } } virtual void onCavity( - EntityArray& oldElements, - EntityArray& newEntities) + ma::EntityArray& oldElements, + ma::EntityArray& newEntities) { others.onCavity(oldElements,newEntities); } @@ -109,7 +108,7 @@ static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f) ma::SolutionTransfer* setBezierSolutionTransfers( const std::vector& fields) { - ma::SolutionTransfer* st = new ma::SolutionTransfers(); + ma::SolutionTransfers* st = new ma::SolutionTransfers(); for (std::size_t i = 0; i < fields.size(); i++) { st->add(createBezierSolutionTransfer(fields[i])); } From 9804ab7d1dab0e41e90241ee46b1606b476ff7a6 Mon Sep 17 00:00:00 2001 From: Avinash Date: Fri, 10 Apr 2020 13:22:35 -0400 Subject: [PATCH 11/36] Bezier representation of field values after Refine. --- crv/crvBezierSolutionTransfer.cc | 31 +++++++++++++++++-------------- ma/maRefine.cc | 8 +++++++- 2 files changed, 24 insertions(+), 15 deletions(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index a141c1878..42a776e71 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -22,6 +22,7 @@ #include #include +#include namespace crv { class CrvBezierSolutionTransfer : public ma::SolutionTransfer @@ -49,22 +50,24 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer ma::EntityArray& newEntities) { others.onRefine(parent,newEntities); - apf::FieldShape *shape = others.shape; apf::Mesh *m = others.mesh; - int type = m->getType(newEntities[0]); - int td = apf::Mesh::typeDimension[type]; - int order = shape->getOrder(); - int n = shape->getEntityShape( - apf::Mesh::simplexTypes[td])->countNodes(); - int ne = shape->countNodesOn( - apf::Mesh::simplexTypes[td]); - apf::NewArray c; - crv::getBezierTransformationCoefficients(order, - apf::Mesh::simplexTypes[td], c); + for (size_t i = 0; i getType(newEntities[i]); + apf::FieldShape *shape = others.shape; + int order = shape->getOrder(); + if (type != 0 && getNumInternalControlPoints(type, order) > 0) { + int td = apf::Mesh::typeDimension[type]; + int n = shape->getEntityShape( + apf::Mesh::simplexTypes[td])->countNodes(); + int ne = shape->countNodesOn( + apf::Mesh::simplexTypes[td]); + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[td], c); - for( size_t i = 0; i < newEntities.getSize(); i++) { - crv::convertInterpolationFieldPoints(newEntities[i], - others.field, n, ne, c); + crv::convertInterpolationFieldPoints(newEntities[i], + others.field, n, ne, c); + } } } virtual void onCavity( diff --git a/ma/maRefine.cc b/ma/maRefine.cc index 4325186fb..c06804d38 100644 --- a/ma/maRefine.cc +++ b/ma/maRefine.cc @@ -345,14 +345,20 @@ void transferElements(Refine* r) Adapt* a = r->adapt; Mesh* m = a->mesh; SolutionTransfer* st = a->solutionTransfer; - int td = st->getTransferDimension(); + int td = a->shape->getTransferDimension(); + for (int d = td; d <= m->getDimension(); ++d) + for (size_t i=0; i < r->toSplit[d].getSize(); ++i) + a->shape->onRefine(r->toSplit[d][i],r->newEntities[d][i]); + td = st->getTransferDimension(); for (int d = td; d <= m->getDimension(); ++d) for (size_t i=0; i < r->toSplit[d].getSize(); ++i) st->onRefine(r->toSplit[d][i],r->newEntities[d][i]); + /* td = a->shape->getTransferDimension(); for (int d = td; d <= m->getDimension(); ++d) for (size_t i=0; i < r->toSplit[d].getSize(); ++i) a->shape->onRefine(r->toSplit[d][i],r->newEntities[d][i]); + */ } void forgetNewEntities(Refine* r) From e10035f168c5844f94952b41f66c70f66f25a5c3 Mon Sep 17 00:00:00 2001 From: Avinash Date: Mon, 13 Apr 2020 11:34:14 -0400 Subject: [PATCH 12/36] Solution transfer: Bezier Field representation on Cavity. --- crv/crvBezierSolutionTransfer.cc | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 42a776e71..e2d8ec39c 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -75,6 +75,25 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer ma::EntityArray& newEntities) { others.onCavity(oldElements,newEntities); + apf::Mesh *m = others.mesh; + for (size_t i = 0; i getType(newEntities[i]); + apf::FieldShape *shape = others.shape; + int order = shape->getOrder(); + if (type != 0 && getNumInternalControlPoints(type, order) > 0) { + int td = apf::Mesh::typeDimension[type]; + int n = shape->getEntityShape( + apf::Mesh::simplexTypes[td])->countNodes(); + int ne = shape->countNodesOn( + apf::Mesh::simplexTypes[td]); + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[td], c); + + crv::convertInterpolationFieldPoints(newEntities[i], + others.field, n, ne, c); + } + } } /* protected: From 7b571c05b27789871463347b086daabbb9a7f1c3 Mon Sep 17 00:00:00 2001 From: Avinash Date: Sun, 19 Apr 2020 13:02:37 -0400 Subject: [PATCH 13/36] solution field transfer upon refine. --- crv/crvBezierSolutionTransfer.cc | 95 +++++++++++++++++++++++++++++++- 1 file changed, 94 insertions(+), 1 deletion(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index e2d8ec39c..a5ac88b89 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -33,6 +33,22 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer CrvBezierSolutionTransfer(apf::Field* f): verts(f),others(f) { + mesh = others.mesh; + shape = others.shape; + int P = mesh->getShape()->getOrder(); + for (int d = 1; d <= mesh->getDimension(); d++) { + int type = apf::Mesh::simplexTypes[d]; + if (!mesh->getShape()->hasNodesIn(d)) + continue; + int n = mesh->getShape()->getEntityShape(type)->countNodes(); + mth::Matrix A(n,n); + Ai[d].resize(n,n); + getBezierTransformationMatrix(type, P, + A, elem_vert_xi[type]); + invertMatrixWithPLU(getNumControlPoints(type,P), A, Ai[d]); + //crv::getBezierTransformationCoeddicients(P,type,coeffs[d]); + //crv::getInternalBezierTransformationCoefficients(mesh, P, 1, + // apf::Mesh::simplexTypes[d],internalCoeffs[d]); } virtual bool hasNodesOn(int dimension) { @@ -45,12 +61,81 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer { verts.onVertex(parent,xi,vert); } + + void getVertParams(int ptype, apf::MeshEntity** parentVerts, + apf::NewArray& midEdgeVerts, apf::MeshEntity* e, + apf::Vector3 params[4]) + { + int npv = apf::Mesh::adjacentCount[ptype][0]; + int ne = apf::Mesh::adjacentCount[ptype][1]; + apf::Downward verts; + + int nv = mesh->getDownward(e,0,verts); + // first check verts + for (int v = 0; v < nv; ++v){ + bool vert = false; + for (int i = 0; i < npv; ++i){ + if(verts[v] == parentVerts[i]){ + params[v] = elem_vert_xi[ptype][i]; + vert = true; + break; + } + } + + if(!vert){ + for (int i = 0; i < ne; ++i){ + if( verts[v] == midEdgeVerts[i] ){ + params[v] = elem_edge_xi[ptype][i]; + break; + } + } + } + } + } + virtual void onRefine( ma::Entity* parent, ma::EntityArray& newEntities) { others.onRefine(parent,newEntities); - apf::Mesh *m = others.mesh; + + int P = shape->getOrder(); + int parentType = mesh->getType(parent); + apf::Downward parentVerts, parentEdges; + + mesh->getDownward(parent, 0, parentVerts); + mesh->getDownward(parent, 1, parentEdges); + + int ne = apf::Mesh::adjacentCount[parentType][1]; + + apf::NewArray midEdgeVerts(ne); + for (int i = 0; i < ne; ++i){ + if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) + midEdgeVerts[i] = ma::findSplitVert(refine,parentEdges[i]); + else + midEdgeVerts[i] = 0; + } + + apf::Element* elem = apf::createElement(f, parent); + apf::NewArray nodes; + + //TODO scalar field values + apf::getVectorNodes(elem,nodes); + + int n = getNumControlPoints(childType,P); + apf::Vector3 vp[4]; + getVertParams(parentType,parentVerts,midEdgeVerts,newEntities[i], vp); + + mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); + for (int j = 0; j < ni; ++j){ + apf::Vector3 value(0,0,0); + for (int k = 0; k < np; ++k) + value += nodes[k]*B(j+n-ni,k); + mesh->setVector(newEntities[i],j,value); + } + + /* + Mesh *m = others.mesh; for (size_t i = 0; i getType(newEntities[i]); apf::FieldShape *shape = others.shape; @@ -69,6 +154,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer others.field, n, ne, c); } } + */ } virtual void onCavity( ma::EntityArray& oldElements, @@ -95,6 +181,13 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer } } } + protected: + ma::Mesh* mesh; + apf::FieldShape* shape; + mth::Matrix Ai[4]; + public: + //apf::NewArray coeffs[4]; + //apf::NewArray internalCoeffs[4]; /* protected: void convertToBezierFields( From 81f419a8f35a0fd4f3bc328da2d663abb146eaad Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Sun, 19 Apr 2020 13:40:22 -0400 Subject: [PATCH 14/36] Adds bezier solution transfers after Adapt is made This is because inside crvBezierSoluitonTransfer we need access to adapt and refine objects. --- crv/crvAdapt.cc | 12 ++++++++++++ crv/crvBezierSolutionTransfer.cc | 5 +++-- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index 3b6baec15..beb10486e 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -194,6 +194,10 @@ static void flagCleaner(crv::Adapt* a) } } +void getAllBezierFields(ma::Mesh* m, std::vector& fields) +{ +} + void adapt(ma::Input* in) { std::string name = in->mesh->getShape()->getName(); @@ -205,6 +209,14 @@ void adapt(ma::Input* in) double t0 = PCU_Time(); ma::validateInput(in); Adapt* a = new Adapt(in); + + // Setting up bezier field transfer for all fields with Bezier shapes + // This is not the cleanest way of doing this! + std::vector allFields; + getAllBezierFields(a->mesh, allFields); + in->solutionTransfer = crv::setBezierSolutionTransfer(allFields, a); + a->solutionTransfer = in->solutionTransfer; + ma::preBalance(a); fixInvalidElements(a); diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index a5ac88b89..976af9288 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -28,10 +28,11 @@ namespace crv { class CrvBezierSolutionTransfer : public ma::SolutionTransfer { public: + Adapt* adapt; ma::LinearTransfer verts; ma::CavityTransfer others; - CrvBezierSolutionTransfer(apf::Field* f): - verts(f),others(f) + CrvBezierSolutionTransfer(apf::Field* f, Adapt* a): + adapt(a),verts(f),others(f) { mesh = others.mesh; shape = others.shape; From 0d4cf81e253c571419bb37c38222cc883ffe28ab Mon Sep 17 00:00:00 2001 From: Avinash Date: Mon, 20 Apr 2020 09:53:51 -0400 Subject: [PATCH 15/36] Debugging crvBezierSolutionFieldTransfer. --- crv/crv.h | 4 +- crv/crvAdapt.cc | 9 ++- crv/crvBezierSolutionTransfer.cc | 103 ++++++++++++------------------- 3 files changed, 50 insertions(+), 66 deletions(-) diff --git a/crv/crv.h b/crv/crv.h index e8f0d1fe5..a7d9c06e5 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -212,8 +212,8 @@ int getTriNodeIndex(int P, int i, int j); int getTetNodeIndex(int P, int i, int j, int k); /** \brief adds bezier solution transfers */ -ma::SolutionTransfer* setBezierSolutionTransfers( - const std::vector& fields); +//ma::SolutionTransfer* setBezierSolutionTransfers( +// const std::vector& fields, crv::Adapt* a); /** \brief crv fail function */ void fail(const char* why) __attribute__((noreturn)); diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index beb10486e..eb88ac3f6 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -7,6 +7,7 @@ #include "crvAdapt.h" #include "crvShape.h" +#include "crvBezierSolutionTransfer.h" #include #include #include @@ -196,6 +197,12 @@ static void flagCleaner(crv::Adapt* a) void getAllBezierFields(ma::Mesh* m, std::vector& fields) { + for (int i = 0; i < m->countFields(); i++) { + apf::FieldShape* fs = apf::getShape(m->getField(i)); + std::string name = fs->getName(); + if (name == std::string("Bezier")) + fields.push_back(m->getField(i)); + } } void adapt(ma::Input* in) @@ -214,7 +221,7 @@ void adapt(ma::Input* in) // This is not the cleanest way of doing this! std::vector allFields; getAllBezierFields(a->mesh, allFields); - in->solutionTransfer = crv::setBezierSolutionTransfer(allFields, a); + in->solutionTransfer = crv::setBezierSolutionTransfers(allFields, a); a->solutionTransfer = in->solutionTransfer; ma::preBalance(a); diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 976af9288..67e5e678f 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -31,11 +31,13 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer Adapt* adapt; ma::LinearTransfer verts; ma::CavityTransfer others; - CrvBezierSolutionTransfer(apf::Field* f, Adapt* a): + CrvBezierSolutionTransfer(apf::Field* field, Adapt* a): adapt(a),verts(f),others(f) { - mesh = others.mesh; - shape = others.shape; + f = field; + refine = adapt->refine; + mesh = adapt->mesh; + shape = apf::getShape(f); int P = mesh->getShape()->getOrder(); for (int d = 1; d <= mesh->getDimension(); d++) { int type = apf::Mesh::simplexTypes[d]; @@ -50,6 +52,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer //crv::getBezierTransformationCoeddicients(P,type,coeffs[d]); //crv::getInternalBezierTransformationCoefficients(mesh, P, 1, // apf::Mesh::simplexTypes[d],internalCoeffs[d]); + } } virtual bool hasNodesOn(int dimension) { @@ -117,46 +120,41 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer midEdgeVerts[i] = 0; } - apf::Element* elem = apf::createElement(f, parent); - apf::NewArray nodes; + int np = shape->getEntityShape(parentType)->countNodes(); - //TODO scalar field values - apf::getVectorNodes(elem,nodes); + apf::Element* elem = apf::createElement(f, parent); + apf::NewArray nodes; - int n = getNumControlPoints(childType,P); - apf::Vector3 vp[4]; - getVertParams(parentType,parentVerts,midEdgeVerts,newEntities[i], vp); + //TODO scalar field values + apf::getVectorNodes(elem,nodes); - mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); - for (int j = 0; j < ni; ++j){ - apf::Vector3 value(0,0,0); - for (int k = 0; k < np; ++k) - value += nodes[k]*B(j+n-ni,k); - mesh->setVector(newEntities[i],j,value); - } + for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++) { + if (!shape->hasNodesIn(d)) continue; + for (size_t i = 0; i < newEntities.getSize(); i++) { + int childType = mesh->getType(newEntities[i]); + if (d != apf::Mesh::typeDimension[childType]) + continue; - /* - Mesh *m = others.mesh; - for (size_t i = 0; i getType(newEntities[i]); - apf::FieldShape *shape = others.shape; - int order = shape->getOrder(); - if (type != 0 && getNumInternalControlPoints(type, order) > 0) { - int td = apf::Mesh::typeDimension[type]; - int n = shape->getEntityShape( - apf::Mesh::simplexTypes[td])->countNodes(); - int ne = shape->countNodesOn( - apf::Mesh::simplexTypes[td]); - apf::NewArray c; - crv::getBezierTransformationCoefficients(order, - apf::Mesh::simplexTypes[td], c); + int ni = shape->countNodesOn(childType); + int n = getNumControlPoints(childType,P); + apf::Vector3 vp[4]; + getVertParams(parentType,parentVerts, + midEdgeVerts,newEntities[i], vp); - crv::convertInterpolationFieldPoints(newEntities[i], - others.field, n, ne, c); + mth::Matrix A(n,np),B(n,n); + getBezierTransformationMatrix(parentType, childType, P,A, vp); + mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); + + for (int j = 0; j < ni; ++j){ + apf::Vector3 value(0,0,0); + for (int k = 0; k < np; ++k) + value += nodes[k]*B(j+n-ni,k); + apf::setVector(f,newEntities[i],j,value); + } } } - */ } + virtual void onCavity( ma::EntityArray& oldElements, ma::EntityArray& newEntities) @@ -183,50 +181,29 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer } } protected: + apf::Field* f; ma::Mesh* mesh; + ma::Refine* refine; apf::FieldShape* shape; mth::Matrix Ai[4]; - public: + //public: //apf::NewArray coeffs[4]; //apf::NewArray internalCoeffs[4]; /* - protected: - void convertToBezierFields( - int dimension, - EntityArray& newEntities) - { - //apf::FieldShape *fs = others.shape; - //std::string name = fs->getName(); - //if (name == std::string("Bezier")) { - int td = apf::Mesh::typeDimension[] - int order = shape->getOrder(); - int n = shape->getEntityShape( - apf::Mesh::simplexTypes[dimension])->countNodes(); - int ne = shape->countNodesOn( - apf::Mesh::simplexTypes[dimension]); - apf::NewArray c; - crv::getBezierTransformationCoefficients(order, - apf::Mesh::simplexTypes[dimension], c); - - for( size_t i = 0; i < newEntities.getSize(); i++) { - crv::convertInterpolationFieldPoints(newEntities[i], - field, n, ne, c); - } - } */ }; -static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f) +static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* a) { - return new CrvBezierSolutionTransfer(f); + return new CrvBezierSolutionTransfer(f, a); } ma::SolutionTransfer* setBezierSolutionTransfers( - const std::vector& fields) + const std::vector& fields, Adapt* a) { ma::SolutionTransfers* st = new ma::SolutionTransfers(); for (std::size_t i = 0; i < fields.size(); i++) { - st->add(createBezierSolutionTransfer(fields[i])); + st->add(createBezierSolutionTransfer(fields[i], a)); } return st; } From 50dfbcc915e15ff48f3be62895e544918de80b09 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Mon, 20 Apr 2020 10:24:40 -0400 Subject: [PATCH 16/36] Forward declares the crv::Adapt class in crv.h --- crv/crv.h | 7 +++++-- crv/crvAdapt.cc | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/crv/crv.h b/crv/crv.h index a7d9c06e5..12a8d4850 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -23,6 +23,9 @@ * \brief the curving functions are contained in this namespace */ namespace crv { +// forward declaration of the crv::Adapt +class Adapt; + /** \brief actually 1 greater than max order */ static unsigned const MAX_ORDER = 19; @@ -212,8 +215,8 @@ int getTriNodeIndex(int P, int i, int j); int getTetNodeIndex(int P, int i, int j, int k); /** \brief adds bezier solution transfers */ -//ma::SolutionTransfer* setBezierSolutionTransfers( -// const std::vector& fields, crv::Adapt* a); +ma::SolutionTransfer* setBezierSolutionTransfers( + const std::vector& fields, crv::Adapt* a); /** \brief crv fail function */ void fail(const char* why) __attribute__((noreturn)); diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index eb88ac3f6..981808208 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -7,7 +7,7 @@ #include "crvAdapt.h" #include "crvShape.h" -#include "crvBezierSolutionTransfer.h" +/* #include "crvBezierSolutionTransfer.h" */ #include #include #include From 06b74c70d718c156492c687a85500e763a541705 Mon Sep 17 00:00:00 2001 From: "Morteza H. Siboni" Date: Mon, 20 Apr 2020 10:25:17 -0400 Subject: [PATCH 17/36] Removes an unnecessary header file --- crv/crvBezierSolutionTransfer.h | 43 --------------------------------- 1 file changed, 43 deletions(-) delete mode 100644 crv/crvBezierSolutionTransfer.h diff --git a/crv/crvBezierSolutionTransfer.h b/crv/crvBezierSolutionTransfer.h deleted file mode 100644 index 9bfd11d75..000000000 --- a/crv/crvBezierSolutionTransfer.h +++ /dev/null @@ -1,43 +0,0 @@ -/* - * Copyright 2015 Scientific Computation Research Center - * - * This work is open source software, licensed under the terms of the - * BSD license as described in the LICENSE file in the top-level directory. - */ - -#ifndef CRVBEZIERSOLUTIONTRANSFER_H -#define CRVBEZIERSOLUTIONTRANSFER_H - -/** \file crvBezierShapes.h - * \brief main file for bezier shape functions */ - -namespace crv { -class CrvBezierSolutionTransfer -{ - public: - virtual ~CrvBezierSolutionTransfer(); - - virtual bool hasNodesOn(int dimension) = 0; - - virtual void onVertex( - apf::MeshElement* parent, - Vector const& xi, - Entity* vert); - - virtual void onRefine( - Entity* parent, - EntityArray& newEntities); - - virtual void onCavity( - EntityArray& oldElements, - EntityArray& newEntities); - - virtual void convertToBezierFields( - int dimension, - EntityArray &newEntities); -} - -CrvBezierSolutionTransfer* createFieldTransfer(f); - - -#endif From dfe6876d358454a776816f154a86c478a6e2dd8e Mon Sep 17 00:00:00 2001 From: Avinash Date: Thu, 7 May 2020 17:11:00 -0400 Subject: [PATCH 18/36] removes commented header file. --- crv/crvAdapt.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index 981808208..330952fd6 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -7,7 +7,6 @@ #include "crvAdapt.h" #include "crvShape.h" -/* #include "crvBezierSolutionTransfer.h" */ #include #include #include From 2c083b749edbc8cb8b53639c68e047370bf05760 Mon Sep 17 00:00:00 2001 From: Avinash Date: Thu, 7 May 2020 20:55:30 -0400 Subject: [PATCH 19/36] fixed solution transfer onCavity and onRefine. --- crv/crvBezierSolutionTransfer.cc | 246 ++++++++++++++++++++++--------- 1 file changed, 177 insertions(+), 69 deletions(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 67e5e678f..4cba54c43 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -25,6 +26,85 @@ #include namespace crv { +static void getVertParams(apf::Mesh* mesh, int ptype, + apf::MeshEntity** parentVerts, + apf::NewArray& midEdgeVerts, + apf::MeshEntity* e, + apf::Vector3 params[4]) +{ + int npv = apf::Mesh::adjacentCount[ptype][0]; + int ne = apf::Mesh::adjacentCount[ptype][1]; + apf::Downward verts; + + int nv = mesh->getDownward(e,0,verts); + // first check verts + for (int v = 0; v < nv; ++v){ + bool vert = false; + for (int i = 0; i < npv; ++i){ + if(verts[v] == parentVerts[i]){ + params[v] = elem_vert_xi[ptype][i]; + vert = true; + break; + } + } + + if(!vert){ + for (int i = 0; i < ne; ++i){ + if( verts[v] == midEdgeVerts[i] ){ + params[v] = elem_edge_xi[ptype][i]; + break; + } + } + } + } +} + +static int getBestElement(int n, + apf::Mesh* mesh, + apf::MeshElement **elems, + apf::Vector3 point, apf::Vector3 &xi) +{ + int iter = 50; + double tol = 1e-4; + + int elemNum = 0; + double value, value2, bestValue = -DBL_MAX; + + apf::Vector3 initialGuess = apf::Vector3(1./4., 1./4., 1./4.); + apf::Vector3 xinew, xyz; + apf::Matrix3x3 Jinv; + + for (int i = 0; i < n; i++) { + apf::Vector3 xin = initialGuess; + + for (int j = 0; j < iter; j++) { + apf::getJacobianInv(elems[i], xin, Jinv); + mapLocalToGlobal(elems[i], xin, xyz); + apf::Vector3 fval = (xyz-point); + + xinew = xin - Jinv*fval; + + if (j > 0 && fval.getLength() > value) break; + + if ( (xinew-xin).getLength() < tol ) { + value2 = ma::getInsideness(mesh, apf::getMeshEntity(elems[i]), xinew); + if ( value2 > bestValue) { + bestValue = value2; + elemNum = i; + xi = xinew; + break; + } + } + else { + value = fval.getLength(); + xin = xinew; + } + } + } + return (elemNum); +} + + class CrvBezierSolutionTransfer : public ma::SolutionTransfer { public: @@ -32,7 +112,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer ma::LinearTransfer verts; ma::CavityTransfer others; CrvBezierSolutionTransfer(apf::Field* field, Adapt* a): - adapt(a),verts(f),others(f) + adapt(a),verts(field),others(field) { f = field; refine = adapt->refine; @@ -49,15 +129,13 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer getBezierTransformationMatrix(type, P, A, elem_vert_xi[type]); invertMatrixWithPLU(getNumControlPoints(type,P), A, Ai[d]); - //crv::getBezierTransformationCoeddicients(P,type,coeffs[d]); - //crv::getInternalBezierTransformationCoefficients(mesh, P, 1, - // apf::Mesh::simplexTypes[d],internalCoeffs[d]); } } virtual bool hasNodesOn(int dimension) { return others.hasNodesOn(dimension); } + virtual void onVertex( apf::MeshElement* parent, ma::Vector const& xi, @@ -66,43 +144,10 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer verts.onVertex(parent,xi,vert); } - void getVertParams(int ptype, apf::MeshEntity** parentVerts, - apf::NewArray& midEdgeVerts, apf::MeshEntity* e, - apf::Vector3 params[4]) - { - int npv = apf::Mesh::adjacentCount[ptype][0]; - int ne = apf::Mesh::adjacentCount[ptype][1]; - apf::Downward verts; - - int nv = mesh->getDownward(e,0,verts); - // first check verts - for (int v = 0; v < nv; ++v){ - bool vert = false; - for (int i = 0; i < npv; ++i){ - if(verts[v] == parentVerts[i]){ - params[v] = elem_vert_xi[ptype][i]; - vert = true; - break; - } - } - - if(!vert){ - for (int i = 0; i < ne; ++i){ - if( verts[v] == midEdgeVerts[i] ){ - params[v] = elem_edge_xi[ptype][i]; - break; - } - } - } - } - } - virtual void onRefine( ma::Entity* parent, ma::EntityArray& newEntities) { - others.onRefine(parent,newEntities); - int P = shape->getOrder(); int parentType = mesh->getType(parent); apf::Downward parentVerts, parentEdges; @@ -114,8 +159,16 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer apf::NewArray midEdgeVerts(ne); for (int i = 0; i < ne; ++i){ - if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) + if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) { midEdgeVerts[i] = ma::findSplitVert(refine,parentEdges[i]); + apf::Element* elemP = apf::createElement(f, parentEdges[i]); + apf::Vector3 xiMid = apf::Vector3(0,0,0); + apf::NewArray val; + val.allocate(apf::countComponents(f)); + apf::getComponents(elemP, xiMid, &(val[0])); + apf::setComponents(f, midEdgeVerts[i], 0, &(val[0])); + apf::destroyElement(elemP); + } else midEdgeVerts[i] = 0; } @@ -123,10 +176,15 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer int np = shape->getEntityShape(parentType)->countNodes(); apf::Element* elem = apf::createElement(f, parent); - apf::NewArray nodes; + apf::NewArray Vnodes; + apf::NewArray Snodes; - //TODO scalar field values - apf::getVectorNodes(elem,nodes); + if (apf::getValueType(f) == apf::VECTOR) { + apf::getVectorNodes(elem,Vnodes); + } + else if (apf::getValueType(f) == apf::SCALAR) { + apf::getScalarNodes(elem,Snodes); + } for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++) { if (!shape->hasNodesIn(d)) continue; @@ -138,7 +196,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer int ni = shape->countNodesOn(childType); int n = getNumControlPoints(childType,P); apf::Vector3 vp[4]; - getVertParams(parentType,parentVerts, + getVertParams(mesh, parentType,parentVerts, midEdgeVerts,newEntities[i], vp); mth::Matrix A(n,np),B(n,n); @@ -146,51 +204,101 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); for (int j = 0; j < ni; ++j){ - apf::Vector3 value(0,0,0); - for (int k = 0; k < np; ++k) - value += nodes[k]*B(j+n-ni,k); - apf::setVector(f,newEntities[i],j,value); + + if (apf::getValueType(f) == apf::VECTOR) { + apf::Vector3 Vvalue(0,0,0); + for (int k = 0; k < np; ++k) + Vvalue += Vnodes[k]*B(j+n-ni,k); + apf::setVector(f,newEntities[i],j,Vvalue); + } + else if (apf::getValueType(f) == apf::SCALAR) { + double Svalue = 0.; + + for (int k = 0; k < np; ++k) + Svalue += Snodes[k]*B(j+n-ni,k); + apf::setScalar(f,newEntities[i],j,Svalue); + } } } } + apf::destroyElement(elem); + } + + void setInterpolatingFieldValues(apf::Mesh* mesh, + ma::EntityArray &oldElements, + apf::MeshEntity* newEnt, int nodeNum, + apf::Vector3 xyz, apf::Field* field) + { + apf::Vector3 xi; + int entNum = -1; + apf::NewArray elems(oldElements.getSize()); + for (size_t i = 0; i < oldElements.getSize(); i++) + elems[i] = apf::createMeshElement(mesh, oldElements[i]); + // For a point outside the cavity getBestElement routine + // returns an extrapolated xi + // should perform a check for negative xi values + entNum = getBestElement(oldElements.getSize(), mesh, + &elems[0], xyz, xi); + + PCU_ALWAYS_ASSERT(entNum >= 0); + + apf::Element* elemP = apf::createElement(f, oldElements[entNum]); + apf::NewArray val; + val.allocate(apf::countComponents(field)); + apf::getComponents(elemP, xi, &(val[0])); + apf::setComponents(f, newEnt, nodeNum, &(val[0])); + + for (size_t i = 0; i < oldElements.getSize(); i++) + apf::destroyMeshElement(elems[i]); + + apf::destroyElement(elemP); } virtual void onCavity( ma::EntityArray& oldElements, ma::EntityArray& newEntities) { - others.onCavity(oldElements,newEntities); - apf::Mesh *m = others.mesh; - for (size_t i = 0; i getType(newEntities[i]); - apf::FieldShape *shape = others.shape; - int order = shape->getOrder(); - if (type != 0 && getNumInternalControlPoints(type, order) > 0) { + + if (getDimension(mesh, oldElements[0]) < 3) + return; + + for (int d = 1; d < 4; d++) { + for(size_t i = 0; i getType(newEntities[i]); int td = apf::Mesh::typeDimension[type]; - int n = shape->getEntityShape( - apf::Mesh::simplexTypes[td])->countNodes(); - int ne = shape->countNodesOn( - apf::Mesh::simplexTypes[td]); - apf::NewArray c; - crv::getBezierTransformationCoefficients(order, - apf::Mesh::simplexTypes[td], c); - - crv::convertInterpolationFieldPoints(newEntities[i], - others.field, n, ne, c); - } + if (td != d) continue; + + int order = shape->getOrder(); + int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); + + if (!ne) continue; + + apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); + apf::Vector3 point, xinew; + int n = shape->getEntityShape(type)->countNodes(); + + for (int j = 0; j < ne; j++) { + shape->getNodeXi(type, j, xinew); + apf::mapLocalToGlobal(me, xinew, point); + setInterpolatingFieldValues(mesh, oldElements, + newEntities[i], j, point, f); + } + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[td], c); + + crv::convertInterpolationFieldPoints(newEntities[i], f, n, ne, c); + apf::destroyMeshElement(me); + } } } + protected: apf::Field* f; ma::Mesh* mesh; ma::Refine* refine; apf::FieldShape* shape; mth::Matrix Ai[4]; - //public: - //apf::NewArray coeffs[4]; - //apf::NewArray internalCoeffs[4]; - /* - */ }; static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* a) From fb2b3ed9ffafd3d73b177f920d0d2ff9e26496b9 Mon Sep 17 00:00:00 2001 From: Avinash Date: Fri, 15 May 2020 18:23:19 -0400 Subject: [PATCH 20/36] Appends the crv Solution transfers to the existing ones. --- crv/crvBezierSolutionTransfer.cc | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 4cba54c43..72005ad0d 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -148,6 +148,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer ma::Entity* parent, ma::EntityArray& newEntities) { + std::cout<<" Inside On Redine of crvSolutionTransfer"<getOrder(); int parentType = mesh->getType(parent); apf::Downward parentVerts, parentEdges; @@ -173,6 +174,7 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer midEdgeVerts[i] = 0; } + std::cout<<" Inside On Refine -- done new nodes due refine "<getEntityShape(parentType)->countNodes(); apf::Element* elem = apf::createElement(f, parent); @@ -309,7 +311,10 @@ static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* ma::SolutionTransfer* setBezierSolutionTransfers( const std::vector& fields, Adapt* a) { - ma::SolutionTransfers* st = new ma::SolutionTransfers(); + ma::SolutionTransfers *st = dynamic_cast(a->solutionTransfer); + if (!st) + st = new ma::SolutionTransfers(); + for (std::size_t i = 0; i < fields.size(); i++) { st->add(createBezierSolutionTransfer(fields[i], a)); } From 613983dc050fc66cf3bda6c1114b4513aa31c1a2 Mon Sep 17 00:00:00 2001 From: Avinash Date: Sat, 23 May 2020 00:17:29 -0400 Subject: [PATCH 21/36] removed print statements. --- ma/maRefine.cc | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/ma/maRefine.cc b/ma/maRefine.cc index c06804d38..08f69c208 100644 --- a/ma/maRefine.cc +++ b/ma/maRefine.cc @@ -20,7 +20,7 @@ #include "maLayer.h" #include #include - +#include namespace ma { void addEdgePreAllocation(Refine* r, Entity* e, int counts[4]) @@ -347,18 +347,16 @@ void transferElements(Refine* r) SolutionTransfer* st = a->solutionTransfer; int td = a->shape->getTransferDimension(); for (int d = td; d <= m->getDimension(); ++d) - for (size_t i=0; i < r->toSplit[d].getSize(); ++i) + for (size_t i=0; i < r->toSplit[d].getSize(); ++i) { + std::cout<<" doing mesh refinement for dimension "<shape->onRefine(r->toSplit[d][i],r->newEntities[d][i]); + } td = st->getTransferDimension(); for (int d = td; d <= m->getDimension(); ++d) - for (size_t i=0; i < r->toSplit[d].getSize(); ++i) + for (size_t i=0; i < r->toSplit[d].getSize(); ++i) { + std::cout<<" doing solution transfer for dimension "<onRefine(r->toSplit[d][i],r->newEntities[d][i]); - /* - td = a->shape->getTransferDimension(); - for (int d = td; d <= m->getDimension(); ++d) - for (size_t i=0; i < r->toSplit[d].getSize(); ++i) - a->shape->onRefine(r->toSplit[d][i],r->newEntities[d][i]); - */ + } } void forgetNewEntities(Refine* r) From ffbffcc4a4f8336ec5eaa2fdb537ebb6511b3d13 Mon Sep 17 00:00:00 2001 From: Avinash Date: Sun, 24 May 2020 12:51:35 -0400 Subject: [PATCH 22/36] getFieldTransferName added. --- crv/crvBezierSolutionTransfer.cc | 20 ++++++++--- ma/maRefine.cc | 2 -- ma/maSolutionTransfer.cc | 58 ++++++++++++++++++++++++++++++-- ma/maSolutionTransfer.h | 6 +++- 4 files changed, 76 insertions(+), 10 deletions(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 72005ad0d..08591a86d 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -109,10 +109,12 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer { public: Adapt* adapt; - ma::LinearTransfer verts; + //ma::LinearTransfer verts; ma::CavityTransfer others; + CrvBezierSolutionTransfer(apf::Field* field, Adapt* a): - adapt(a),verts(field),others(field) + adapt(a),others(field) + //adapt(a),verts(field),others(field) { f = field; refine = adapt->refine; @@ -136,6 +138,11 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer return others.hasNodesOn(dimension); } + virtual const char* getTransferFieldName() { + return apf::getName(f); + } + + /* virtual void onVertex( apf::MeshElement* parent, ma::Vector const& xi, @@ -143,12 +150,11 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer { verts.onVertex(parent,xi,vert); } - +*/ virtual void onRefine( ma::Entity* parent, ma::EntityArray& newEntities) { - std::cout<<" Inside On Redine of crvSolutionTransfer"<getOrder(); int parentType = mesh->getType(parent); apf::Downward parentVerts, parentEdges; @@ -174,7 +180,6 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer midEdgeVerts[i] = 0; } - std::cout<<" Inside On Refine -- done new nodes due refine "<getEntityShape(parentType)->countNodes(); apf::Element* elem = apf::createElement(f, parent); @@ -318,6 +323,11 @@ ma::SolutionTransfer* setBezierSolutionTransfers( for (std::size_t i = 0; i < fields.size(); i++) { st->add(createBezierSolutionTransfer(fields[i], a)); } + std::vector trans = st->transfers; + for (std::size_t i = 0; i < trans.size(); i++) { + const char* name = trans[i]->getTransferFieldName(); + std::cout<<" field name added "<< name<shape->getTransferDimension(); for (int d = td; d <= m->getDimension(); ++d) for (size_t i=0; i < r->toSplit[d].getSize(); ++i) { - std::cout<<" doing mesh refinement for dimension "<shape->onRefine(r->toSplit[d][i],r->newEntities[d][i]); } td = st->getTransferDimension(); for (int d = td; d <= m->getDimension(); ++d) for (size_t i=0; i < r->toSplit[d].getSize(); ++i) { - std::cout<<" doing solution transfer for dimension "<onRefine(r->toSplit[d][i],r->newEntities[d][i]); } } diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index e7ae9f403..999ea3144 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -10,7 +10,8 @@ #include "maSolutionTransfer.h" #include "maAffine.h" #include "maMap.h" - +#include "apfField.h" +#include namespace ma { SolutionTransfer::~SolutionTransfer() @@ -60,6 +61,10 @@ int SolutionTransfer::getTransferDimension() return transferDimension; } +const char* SolutionTransfer::getTransferFieldName() +{ + return 0; +} FieldTransfer::FieldTransfer(apf::Field* f) { @@ -74,6 +79,11 @@ bool FieldTransfer::hasNodesOn(int dimension) return shape->hasNodesIn(dimension); } +const char* FieldTransfer::getTransferFieldName() +{ + return apf::getName(field); +} + void LinearTransfer::onVertex( apf::MeshElement* parent, Vector const& xi, @@ -172,7 +182,45 @@ void CavityTransfer::onCavity( EntityArray& oldElements, EntityArray& newEntities) { + apf::FieldShape *fs = apf::getShape(field); + std::string name = fs->getName(); + std::cout<<" field name inside oncavity "<getType(newEntities[i]); + if (type == 4) { + double val = apf::getScalar(field, newEntities[i], 0); + //double val = apf::measure(mesh, newEntities[i]); + sumNew += val; + } + } + for (size_t i = 0; i < newEntities.getSize(); i++) { + int type = mesh->getType(newEntities[i]); + if (type == 4) { + double val = apf::getScalar(field, newEntities[i], 0); + //double val = apf::measure(mesh, newEntities[i]); + //apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); + apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); + } + } + } + std::cout<<" value after adapt "<countFields(); ++i) { apf::Field* f = m->getField(i); - this->add(createFieldTransfer(f)); + std::string name = f->getShape()->getName(); + if (name != std::string("Bezier")) + this->add(createFieldTransfer(f)); } } diff --git a/ma/maSolutionTransfer.h b/ma/maSolutionTransfer.h index 1d59e65b4..5e1f11a83 100644 --- a/ma/maSolutionTransfer.h +++ b/ma/maSolutionTransfer.h @@ -15,6 +15,7 @@ #include #include +#include #include #include #include "maMesh.h" @@ -76,6 +77,7 @@ class SolutionTransfer EntityArray& newEntities); /** \brief for internal MeshAdapt use */ int getTransferDimension(); + virtual const char* getTransferFieldName(); }; class FieldTransfer : public SolutionTransfer @@ -87,6 +89,7 @@ class FieldTransfer : public SolutionTransfer apf::Mesh* mesh; apf::FieldShape* shape; apf::NewArray value; + virtual const char* getTransferFieldName(); }; class LinearTransfer : public FieldTransfer @@ -94,6 +97,7 @@ class LinearTransfer : public FieldTransfer public: LinearTransfer(apf::Field* f): FieldTransfer(f) {} + apf::Field* field; virtual void onVertex( apf::MeshElement* parent, Vector const& xi, @@ -161,7 +165,7 @@ class SolutionTransfers : public SolutionTransfer virtual void onCavity( EntityArray& oldElements, EntityArray& newEntities); - private: + //private typedef std::vector Transfers; Transfers transfers; }; From e273e25ff2725bab174d0d044b5c6543585f2b28 Mon Sep 17 00:00:00 2001 From: Avinash Date: Wed, 27 May 2020 22:29:25 -0400 Subject: [PATCH 23/36] constant field transfer in a cavity, weighted by old to new values. --- ma/maSolutionTransfer.cc | 16 +++++++++++----- ma/maSolutionTransfer.h | 1 - 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index 999ea3144..94ad10fa0 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -184,7 +184,6 @@ void CavityTransfer::onCavity( { apf::FieldShape *fs = apf::getShape(field); std::string name = fs->getName(); - std::cout<<" field name inside oncavity "<getType(newEntities[i]); if (type == 4) { double val = apf::getScalar(field, newEntities[i], 0); //double val = apf::measure(mesh, newEntities[i]); - sumNew += val; + sumNew += val*apf::measure(mesh, newEntities[i]); } } for (size_t i = 0; i < newEntities.getSize(); i++) { @@ -219,8 +218,15 @@ void CavityTransfer::onCavity( apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); } } + for (size_t i = 0; i < newEntities.getSize(); i++) { + int type = mesh->getType(newEntities[i]); + if (type == 4) { + double val = apf::getScalar(field, newEntities[i], 0); + //double val = apf::measure(mesh, newEntities[i]); + sumNewVal += val*apf::measure(mesh, newEntities[i]); + } + } } - std::cout<<" value after adapt "< Date: Sat, 6 Jun 2020 05:42:30 -0400 Subject: [PATCH 24/36] bezier fields + st. --- crv/crvBezierFields.cc | 53 ++- crv/crvBezierSolutionTransfer.cc | 541 ++++++++++++++++--------------- crv/crvShapeHandler.cc | 2 +- ma/maSolutionTransfer.cc | 1 + 4 files changed, 324 insertions(+), 273 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index 06cabcfd6..c63a700b7 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -16,12 +16,12 @@ #include #include #include - +#include #include namespace crv { -static void convertVectorFieldInterpolationPoints(int n, int ne, +static void convertVectorFieldInterpolationPoints(int n, int ne, apf::NewArray& nodes, apf::NewArray& c, apf::NewArray& newNodes){ @@ -34,10 +34,11 @@ static void convertVectorFieldInterpolationPoints(int n, int ne, newNodes[i] += nodes[j]*c[i*n+j]; } -static void convertScalarFieldInterpolationPoints(int n, int ne, +static void convertScalarFieldInterpolationPoints(int n, int ne, apf::NewArray& nodes, apf::NewArray& c, - apf::NewArray& newNodes){ + apf::NewArray& newNodes) +{ for(int i = 0; i < ne; ++i) newNodes[i] = 0.; @@ -49,7 +50,8 @@ static void convertScalarFieldInterpolationPoints(int n, int ne, void convertInterpolationFieldPoints(apf::MeshEntity* e, apf::Field* f, - int n, int ne, apf::NewArray& c){ + int n, int ne, apf::NewArray& c) +{ apf::NewArray l, b(ne); apf::NewArray ls, bs(ne); @@ -82,11 +84,16 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) // Zero out the blending orders (since they are used internally // in Bezier shapes (triangles and tets) but save the old ones to // revert back to at the end of this (so the rest of Bezier related - // procedures work as intended) - //int oldBlendings[apf::Mesh::TYPES]; - //for (int i = 0; i < apf::Mesh::TYPES; i++) - // oldBlendings[i] = getBlendingOrder(i); - setBlendingOrder(apf::Mesh::TYPES, 0); + // procedures work as intendedi) + // + /* + int oldBlendings[apf::Mesh::TYPES]; + for (int i = 0; i < apf::Mesh::TYPES; i++) { + oldBlendings[i] = getBlendingOrder(i); + std::cout<<" blending order "<< oldBlendings[i] <<" for type "<< i<getOrder(); @@ -96,13 +103,14 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) //transferFields(m_mesh, f, fnew); int md = m_mesh->getDimension(); - //int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]); + // + int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]); //blendingOrder = 0; // go downward, and convert interpolating to control points - //int startDim = md - (blendingOrder > 0); + int startDim = md - (blendingOrder > 0); - //for(int d = startDim; d >= 1; --d){ - for(int d = md; d >= 1; d--){ + for(int d = startDim; d >= 1; --d){ + //for(int d = md; d >= 1; d--){ if(!fs->hasNodesIn(d)) continue; int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); @@ -118,6 +126,23 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) m_mesh->end(it); } + for( int d = 2; d <= md; ++d){ + if(!fs->hasNodesIn(d) || + getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; + int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); + int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); + apf::NewArray c; + getInternalBezierTransformationCoefficients(m_mesh,order,1, + apf::Mesh::simplexTypes[d],c); + apf::MeshEntity* e; + apf::MeshIterator* it = m_mesh->begin(d); + while ((e = m_mesh->iterate(it))){ + if(!isBoundaryEntity(m_mesh,e) && m_mesh->isOwned(e)) + convertInterpolationFieldPoints(e,f,n-ne,ne,c); + } + m_mesh->end(it); + } + apf::synchronize(f); } diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 08591a86d..e71c17a9b 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -1,5 +1,5 @@ /* -1;3409;0c * Copyright 2015 Scientific Computation Research Center + 1;3409;0c * Copyright 2015 Scientific Computation Research Center * * This work is open source software, licensed under the terms of the * BSD license as described in the LICENSE file in the top-level directory. @@ -26,95 +26,95 @@ #include namespace crv { -static void getVertParams(apf::Mesh* mesh, int ptype, - apf::MeshEntity** parentVerts, - apf::NewArray& midEdgeVerts, - apf::MeshEntity* e, - apf::Vector3 params[4]) -{ - int npv = apf::Mesh::adjacentCount[ptype][0]; - int ne = apf::Mesh::adjacentCount[ptype][1]; - apf::Downward verts; - - int nv = mesh->getDownward(e,0,verts); - // first check verts - for (int v = 0; v < nv; ++v){ - bool vert = false; - for (int i = 0; i < npv; ++i){ - if(verts[v] == parentVerts[i]){ - params[v] = elem_vert_xi[ptype][i]; - vert = true; - break; + static void getVertParams(apf::Mesh* mesh, int ptype, + apf::MeshEntity** parentVerts, + apf::NewArray& midEdgeVerts, + apf::MeshEntity* e, + apf::Vector3 params[4]) + { + int npv = apf::Mesh::adjacentCount[ptype][0]; + int ne = apf::Mesh::adjacentCount[ptype][1]; + apf::Downward verts; + + int nv = mesh->getDownward(e,0,verts); + // first check verts + for (int v = 0; v < nv; ++v){ + bool vert = false; + for (int i = 0; i < npv; ++i){ + if(verts[v] == parentVerts[i]){ + params[v] = elem_vert_xi[ptype][i]; + vert = true; + break; + } } - } - if(!vert){ - for (int i = 0; i < ne; ++i){ - if( verts[v] == midEdgeVerts[i] ){ - params[v] = elem_edge_xi[ptype][i]; - break; - } + if(!vert){ + for (int i = 0; i < ne; ++i){ + if( verts[v] == midEdgeVerts[i] ){ + params[v] = elem_edge_xi[ptype][i]; + break; + } + } } } } -} - -static int getBestElement(int n, - apf::Mesh* mesh, - apf::MeshElement **elems, - apf::Vector3 point, apf::Vector3 &xi) -{ - int iter = 50; - double tol = 1e-4; - - int elemNum = 0; - double value, value2, bestValue = -DBL_MAX; - - apf::Vector3 initialGuess = apf::Vector3(1./4., 1./4., 1./4.); - apf::Vector3 xinew, xyz; - apf::Matrix3x3 Jinv; - - for (int i = 0; i < n; i++) { - apf::Vector3 xin = initialGuess; - - for (int j = 0; j < iter; j++) { - apf::getJacobianInv(elems[i], xin, Jinv); - mapLocalToGlobal(elems[i], xin, xyz); - apf::Vector3 fval = (xyz-point); - - xinew = xin - Jinv*fval; - - if (j > 0 && fval.getLength() > value) break; - - if ( (xinew-xin).getLength() < tol ) { - value2 = ma::getInsideness(mesh, apf::getMeshEntity(elems[i]), xinew); - if ( value2 > bestValue) { - bestValue = value2; - elemNum = i; - xi = xinew; - break; - } - } - else { - value = fval.getLength(); - xin = xinew; + + static int getBestElement(int n, + apf::Mesh* mesh, + apf::MeshElement **elems, + apf::Vector3 point, apf::Vector3 &xi) + { + int iter = 50; + double tol = 1e-4; + + int elemNum = 0; + double value, value2, bestValue = -DBL_MAX; + + apf::Vector3 initialGuess = apf::Vector3(1./4., 1./4., 1./4.); + apf::Vector3 xinew, xyz; + apf::Matrix3x3 Jinv; + + for (int i = 0; i < n; i++) { + apf::Vector3 xin = initialGuess; + + for (int j = 0; j < iter; j++) { + apf::getJacobianInv(elems[i], xin, Jinv); + mapLocalToGlobal(elems[i], xin, xyz); + apf::Vector3 fval = (xyz-point); + + xinew = xin - Jinv*fval; + + if (j > 0 && fval.getLength() > value) break; + + if ( (xinew-xin).getLength() < tol ) { + value2 = ma::getInsideness(mesh, apf::getMeshEntity(elems[i]), xinew); + if ( value2 > bestValue) { + bestValue = value2; + elemNum = i; + xi = xinew; + break; + } + } + else { + value = fval.getLength(); + xin = xinew; + } } } + return (elemNum); } - return (elemNum); -} -class CrvBezierSolutionTransfer : public ma::SolutionTransfer -{ - public: - Adapt* adapt; - //ma::LinearTransfer verts; - ma::CavityTransfer others; + class CrvBezierSolutionTransfer : public ma::SolutionTransfer + { + public: + Adapt* adapt; + //ma::LinearTransfer verts; + ma::CavityTransfer others; - CrvBezierSolutionTransfer(apf::Field* field, Adapt* a): - adapt(a),others(field) - //adapt(a),verts(field),others(field) + CrvBezierSolutionTransfer(apf::Field* field, Adapt* a): + adapt(a),others(field) + //adapt(a),verts(field),others(field) { f = field; refine = adapt->refine; @@ -122,10 +122,10 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer shape = apf::getShape(f); int P = mesh->getShape()->getOrder(); for (int d = 1; d <= mesh->getDimension(); d++) { - int type = apf::Mesh::simplexTypes[d]; - if (!mesh->getShape()->hasNodesIn(d)) - continue; - int n = mesh->getShape()->getEntityShape(type)->countNodes(); + int type = apf::Mesh::simplexTypes[d]; + if (!mesh->getShape()->hasNodesIn(d)) + continue; + int n = mesh->getShape()->getEntityShape(type)->countNodes(); mth::Matrix A(n,n); Ai[d].resize(n,n); getBezierTransformationMatrix(type, P, @@ -133,202 +133,227 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer invertMatrixWithPLU(getNumControlPoints(type,P), A, Ai[d]); } } - virtual bool hasNodesOn(int dimension) - { - return others.hasNodesOn(dimension); - } + virtual bool hasNodesOn(int dimension) + { + return others.hasNodesOn(dimension); + } - virtual const char* getTransferFieldName() { - return apf::getName(f); - } + virtual const char* getTransferFieldName() { + return apf::getName(f); + } - /* - virtual void onVertex( - apf::MeshElement* parent, - ma::Vector const& xi, - ma::Entity* vert) - { - verts.onVertex(parent,xi,vert); - } -*/ - virtual void onRefine( - ma::Entity* parent, - ma::EntityArray& newEntities) - { - int P = shape->getOrder(); - int parentType = mesh->getType(parent); - apf::Downward parentVerts, parentEdges; - - mesh->getDownward(parent, 0, parentVerts); - mesh->getDownward(parent, 1, parentEdges); - - int ne = apf::Mesh::adjacentCount[parentType][1]; - - apf::NewArray midEdgeVerts(ne); - for (int i = 0; i < ne; ++i){ - if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) { - midEdgeVerts[i] = ma::findSplitVert(refine,parentEdges[i]); - apf::Element* elemP = apf::createElement(f, parentEdges[i]); - apf::Vector3 xiMid = apf::Vector3(0,0,0); - apf::NewArray val; - val.allocate(apf::countComponents(f)); - apf::getComponents(elemP, xiMid, &(val[0])); - apf::setComponents(f, midEdgeVerts[i], 0, &(val[0])); - apf::destroyElement(elemP); - } - else - midEdgeVerts[i] = 0; + virtual void onVertex( + apf::MeshElement* parent, + ma::Vector const& xi, + ma::Entity* vert) + { + apf::Element* e = apf::createElement(f, parent); + apf::NewArray value; + value.allocate(apf::countComponents(f)); + apf::getComponents(e, xi, &(value[0])); + apf::setComponents(f, vert, 0, &(value[0])); + apf::destroyElement(e); } - int np = shape->getEntityShape(parentType)->countNodes(); + virtual void onRefine( + ma::Entity* parent, + ma::EntityArray& newEntities) + { + int P = shape->getOrder(); + int parentType = mesh->getType(parent); + apf::Downward parentVerts, parentEdges; - apf::Element* elem = apf::createElement(f, parent); - apf::NewArray Vnodes; - apf::NewArray Snodes; + mesh->getDownward(parent, 0, parentVerts); + mesh->getDownward(parent, 1, parentEdges); - if (apf::getValueType(f) == apf::VECTOR) { - apf::getVectorNodes(elem,Vnodes); - } - else if (apf::getValueType(f) == apf::SCALAR) { - apf::getScalarNodes(elem,Snodes); - } + int ne = apf::Mesh::adjacentCount[parentType][1]; - for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++) { - if (!shape->hasNodesIn(d)) continue; - for (size_t i = 0; i < newEntities.getSize(); i++) { - int childType = mesh->getType(newEntities[i]); - if (d != apf::Mesh::typeDimension[childType]) - continue; - - int ni = shape->countNodesOn(childType); - int n = getNumControlPoints(childType,P); - apf::Vector3 vp[4]; - getVertParams(mesh, parentType,parentVerts, - midEdgeVerts,newEntities[i], vp); - - mth::Matrix A(n,np),B(n,n); - getBezierTransformationMatrix(parentType, childType, P,A, vp); - mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); - - for (int j = 0; j < ni; ++j){ - - if (apf::getValueType(f) == apf::VECTOR) { - apf::Vector3 Vvalue(0,0,0); - for (int k = 0; k < np; ++k) - Vvalue += Vnodes[k]*B(j+n-ni,k); - apf::setVector(f,newEntities[i],j,Vvalue); - } - else if (apf::getValueType(f) == apf::SCALAR) { - double Svalue = 0.; + //std::cout<<" field name "< midEdgeVerts(ne); + + for (int i = 0; i < ne; ++i){ + if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) + midEdgeVerts[i] = ma::findSplitVert(refine,parentEdges[i]); + else + midEdgeVerts[i] = 0; + } - for (int k = 0; k < np; ++k) - Svalue += Snodes[k]*B(j+n-ni,k); - apf::setScalar(f,newEntities[i],j,Svalue); + int np = shape->getEntityShape(parentType)->countNodes(); + + apf::Element* elem = apf::createElement(f, parent); + apf::NewArray Vnodes; + apf::NewArray Snodes; + + if (apf::getValueType(f) == apf::VECTOR) { + apf::getVectorNodes(elem,Vnodes); + } + else if (apf::getValueType(f) == apf::SCALAR) { + apf::getScalarNodes(elem,Snodes); + } + + for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++) { + if (!shape->hasNodesIn(d)) continue; + for (size_t i = 0; i < newEntities.getSize(); i++) { + int childType = mesh->getType(newEntities[i]); + if (d != apf::Mesh::typeDimension[childType]) + continue; + + //int n = getNumControlPoints(childType,P); + int n = shape->getEntityShape(childType)->countNodes(); + int ni = shape->countNodesOn(childType); + apf::Vector3 vp[4]; + getVertParams(mesh, parentType,parentVerts, + midEdgeVerts,newEntities[i], vp); + + mth::Matrix A(n,np),B(n,np); + getBezierTransformationMatrix(parentType, childType, P,A, vp); + mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); + + for (int j = 0; j < ni; ++j){ + + if (apf::getValueType(f) == apf::VECTOR) { + apf::Vector3 Vvalue(0,0,0); + for (int k = 0; k < np; ++k) + Vvalue += Vnodes[k]*B(j+n-ni,k); + apf::setVector(f,newEntities[i],ni-j-1,Vvalue); + } + else if (apf::getValueType(f) == apf::SCALAR) { + double Svalue = 0.; + + for (int k = 0; k < np; ++k) + Svalue += Snodes[k]*B(j+n-ni,k); + apf::setScalar(f,newEntities[i],ni-j-1,Svalue); + } } - } - } + } + } +/* + //for (int d = apf::Mesh::typeDimension[parentType]; d >= 1; --d){ + for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++){ + if (!shape->hasNodesIn(d)) continue; + + for(size_t i = 0; i getType(newEntities[i]); + int td = apf::Mesh::typeDimension[type]; + if (td != d) continue; + + int ni = shape->countNodesOn(apf::Mesh::simplexTypes[td]); + + if (!ni) continue; + int n = shape->getEntityShape(type)->countNodes(); + + apf::NewArray c; + crv::getBezierTransformationCoefficients(P, + apf::Mesh::simplexTypes[td], c); + crv::convertInterpolationFieldPoints(newEntities[i], + f, n, ni, c); + } + } + */ + apf::destroyElement(elem); } - apf::destroyElement(elem); - } - void setInterpolatingFieldValues(apf::Mesh* mesh, - ma::EntityArray &oldElements, - apf::MeshEntity* newEnt, int nodeNum, - apf::Vector3 xyz, apf::Field* field) - { - apf::Vector3 xi; - int entNum = -1; - apf::NewArray elems(oldElements.getSize()); - for (size_t i = 0; i < oldElements.getSize(); i++) - elems[i] = apf::createMeshElement(mesh, oldElements[i]); - // For a point outside the cavity getBestElement routine - // returns an extrapolated xi - // should perform a check for negative xi values - entNum = getBestElement(oldElements.getSize(), mesh, - &elems[0], xyz, xi); - - PCU_ALWAYS_ASSERT(entNum >= 0); - - apf::Element* elemP = apf::createElement(f, oldElements[entNum]); - apf::NewArray val; - val.allocate(apf::countComponents(field)); - apf::getComponents(elemP, xi, &(val[0])); - apf::setComponents(f, newEnt, nodeNum, &(val[0])); - - for (size_t i = 0; i < oldElements.getSize(); i++) - apf::destroyMeshElement(elems[i]); - - apf::destroyElement(elemP); - } + void setInterpolatingFieldValues(apf::Mesh* mesh, + ma::EntityArray &oldElements, + apf::MeshEntity* newEnt, int nodeNum, + apf::Vector3 xyz, apf::Field* field) + { + apf::Vector3 xi; + int entNum = -1; + int n = oldElements.getSize(); + + apf::NewArray elems(n); + for (int i = 0; i < n; i++) + elems[i] = apf::createMeshElement(mesh, oldElements[i]); + // For a point outside the cavity getBestElement routine + // returns an extrapolated xi + // should perform a check for negative xi values + entNum = getBestElement(n, mesh, + &elems[0], xyz, xi); + + PCU_ALWAYS_ASSERT(entNum >= 0); + + apf::Element* elemP = apf::createElement(f, oldElements[entNum]); + apf::NewArray val; + val.allocate(apf::countComponents(field)); + apf::getComponents(elemP, xi, &(val[0])); + apf::setComponents(f, newEnt, nodeNum, &(val[0])); + + for (int i = 0; i < n; i++) + apf::destroyMeshElement(elems[i]); + + apf::destroyElement(elemP); + } - virtual void onCavity( - ma::EntityArray& oldElements, - ma::EntityArray& newEntities) - { + virtual void onCavity( + ma::EntityArray& oldElements, + ma::EntityArray& newEntities) + { - if (getDimension(mesh, oldElements[0]) < 3) - return; + if (getDimension(mesh, oldElements[0]) < 3) + return; - for (int d = 1; d < 4; d++) { - for(size_t i = 0; i getType(newEntities[i]); - int td = apf::Mesh::typeDimension[type]; - if (td != d) continue; + for (int d = 1; d < 4; d++) { + for(size_t i = 0; i getType(newEntities[i]); + int td = apf::Mesh::typeDimension[type]; + if (td != d) continue; - int order = shape->getOrder(); - int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); + int order = shape->getOrder(); + int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); - if (!ne) continue; + if (!ne) continue; - apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); - apf::Vector3 point, xinew; - int n = shape->getEntityShape(type)->countNodes(); + apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); + apf::Vector3 point, xinew; + int n = shape->getEntityShape(type)->countNodes(); - for (int j = 0; j < ne; j++) { - shape->getNodeXi(type, j, xinew); - apf::mapLocalToGlobal(me, xinew, point); - setInterpolatingFieldValues(mesh, oldElements, - newEntities[i], j, point, f); + for (int j = 0; j < ne; j++) { + shape->getNodeXi(type, j, xinew); + apf::mapLocalToGlobal(me, xinew, point); + setInterpolatingFieldValues(mesh, oldElements, + newEntities[i], j, point, f); + } + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[td], c); + + crv::convertInterpolationFieldPoints(newEntities[i], f, n, ne, c); + apf::destroyMeshElement(me); } - apf::NewArray c; - crv::getBezierTransformationCoefficients(order, - apf::Mesh::simplexTypes[td], c); + } + } - crv::convertInterpolationFieldPoints(newEntities[i], f, n, ne, c); - apf::destroyMeshElement(me); + protected: + apf::Field* f; + ma::Mesh* mesh; + ma::Refine* refine; + apf::FieldShape* shape; + mth::Matrix Ai[4]; + }; + + static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* a) + { + return new CrvBezierSolutionTransfer(f, a); + } + + ma::SolutionTransfer* setBezierSolutionTransfers( + const std::vector& fields, Adapt* a) + { + ma::SolutionTransfers *st = dynamic_cast(a->solutionTransfer); + if (!st) + st = new ma::SolutionTransfers(); + + for (std::size_t i = 0; i < fields.size(); i++) { + st->add(createBezierSolutionTransfer(fields[i], a)); + } + std::vector trans = st->transfers; + for (std::size_t i = 0; i < trans.size(); i++) { + const char* name = trans[i]->getTransferFieldName(); + std::cout<<" field name added "<< name< Ai[4]; -}; - -static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* a) -{ - return new CrvBezierSolutionTransfer(f, a); -} - -ma::SolutionTransfer* setBezierSolutionTransfers( - const std::vector& fields, Adapt* a) -{ - ma::SolutionTransfers *st = dynamic_cast(a->solutionTransfer); - if (!st) - st = new ma::SolutionTransfers(); - - for (std::size_t i = 0; i < fields.size(); i++) { - st->add(createBezierSolutionTransfer(fields[i], a)); - } - std::vector trans = st->transfers; - for (std::size_t i = 0; i < trans.size(); i++) { - const char* name = trans[i]->getTransferFieldName(); - std::cout<<" field name added "<< name<setPoint(newEntities[i],j,point); } } else { - mth::Matrix A(n,np),B(n,n); + mth::Matrix A(n,np),B(n,np); getBezierTransformationMatrix(parentType,childType,P,A,vp); mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); for (int j = 0; j < ni; ++j){ diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index 94ad10fa0..749b774c7 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -215,6 +215,7 @@ void CavityTransfer::onCavity( double val = apf::getScalar(field, newEntities[i], 0); //double val = apf::measure(mesh, newEntities[i]); //apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); + if (val == 0) continue; apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); } } From 44b059c6155997f8662731425f6518eb4485d8f7 Mon Sep 17 00:00:00 2001 From: Avinash Date: Fri, 12 Jun 2020 16:55:01 -0400 Subject: [PATCH 25/36] constant Field cavity transfer. --- ma/CMakeLists.txt | 2 ++ ma/maSolutionTransfer.cc | 18 ++++++++++++++++-- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/ma/CMakeLists.txt b/ma/CMakeLists.txt index 68e72d2ea..6ef1bc036 100644 --- a/ma/CMakeLists.txt +++ b/ma/CMakeLists.txt @@ -56,6 +56,8 @@ set(HEADERS ma.h maInput.h maMesh.h + maMap.h + maAffine.h maSize.h maShape.h maTables.h diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index 749b774c7..bdfcbcea9 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -176,8 +176,21 @@ void CavityTransfer::onRefine( Entity* parent, EntityArray& newEntities) { - transfer(1,&parent,newEntities); + apf::FieldShape *fs = apf::getShape(field); + std::string name = fs->getName(); + if (name != std::string("Constant_3")) + transfer(1,&parent,newEntities); + else { + for (size_t i = 0; i < newEntities.getSize(); i++) { + int type = mesh->getType(newEntities[i]); + if (type == 4) { + double val = apf::getScalar(field, parent, 0); + apf::setScalar(field, newEntities[i], 0, val); + } + } + } } + void CavityTransfer::onCavity( EntityArray& oldElements, EntityArray& newEntities) @@ -193,9 +206,9 @@ void CavityTransfer::onCavity( for (size_t i = 0; i < oldElements.getSize(); i++) { double val = apf::getScalar(field, oldElements[i], 0); sumOld += val*apf::measure(mesh, oldElements[i]); + //sumOld += val; } } - transfer(oldElements.getSize(),&(oldElements[0]),newEntities); double sumNew = 0.; @@ -207,6 +220,7 @@ void CavityTransfer::onCavity( double val = apf::getScalar(field, newEntities[i], 0); //double val = apf::measure(mesh, newEntities[i]); sumNew += val*apf::measure(mesh, newEntities[i]); + //sumNew += val; } } for (size_t i = 0; i < newEntities.getSize(); i++) { From 0b352bfed96f99ef4565191d9f47d9d6fb2c12ad Mon Sep 17 00:00:00 2001 From: Avinash Date: Mon, 15 Jun 2020 10:27:31 -0400 Subject: [PATCH 26/36] index fixed in ST OnRefine. --- crv/crvBezierSolutionTransfer.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index e71c17a9b..3b77fa128 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -216,14 +216,14 @@ namespace crv { apf::Vector3 Vvalue(0,0,0); for (int k = 0; k < np; ++k) Vvalue += Vnodes[k]*B(j+n-ni,k); - apf::setVector(f,newEntities[i],ni-j-1,Vvalue); + apf::setVector(f,newEntities[i],j,Vvalue); } else if (apf::getValueType(f) == apf::SCALAR) { double Svalue = 0.; for (int k = 0; k < np; ++k) Svalue += Snodes[k]*B(j+n-ni,k); - apf::setScalar(f,newEntities[i],ni-j-1,Svalue); + apf::setScalar(f,newEntities[i],j,Svalue); } } } From b77a82132062e6cdf12bccd0472bcfdfaefab84c Mon Sep 17 00:00:00 2001 From: Avinash Date: Fri, 26 Jun 2020 04:06:21 -0400 Subject: [PATCH 27/36] changes specific to Laghos problem. --- crv/crvBezierSolutionTransfer.cc | 121 +++++++++++++++++++++++++++---- crv/crvCurveMesh.cc | 4 +- crv/crvShape.cc | 21 +++--- ma/maShape.cc | 1 + ma/maSolutionTransfer.cc | 2 +- 5 files changed, 123 insertions(+), 26 deletions(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 3b77fa128..3984a01e9 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -24,6 +24,7 @@ #include #include +#include namespace crv { static void getVertParams(apf::Mesh* mesh, int ptype, @@ -64,46 +65,86 @@ namespace crv { apf::MeshElement **elems, apf::Vector3 point, apf::Vector3 &xi) { - int iter = 50; + int iter = 1000; double tol = 1e-4; int elemNum = 0; - double value, value2, bestValue = -DBL_MAX; + double value2, bestValue = -DBL_MAX; apf::Vector3 initialGuess = apf::Vector3(1./4., 1./4., 1./4.); apf::Vector3 xinew, xyz; apf::Matrix3x3 Jinv; + apf::Vector3 fval; + //printf(" inside find best element loop %d \n", n); for (int i = 0; i < n; i++) { apf::Vector3 xin = initialGuess; + apf::Element* e = apf::createElement(mesh->getCoordinateField(), elems[i]); for (int j = 0; j < iter; j++) { apf::getJacobianInv(elems[i], xin, Jinv); - mapLocalToGlobal(elems[i], xin, xyz); - apf::Vector3 fval = (xyz-point); + //mapLocalToGlobal(elems[i], xin, xyz); + apf::getVector(e, xin, xyz); + fval = (xyz-point); - xinew = xin - Jinv*fval; - - if (j > 0 && fval.getLength() > value) break; + xinew = xin - transpose(Jinv)*fval; if ( (xinew-xin).getLength() < tol ) { + //elemNum = i; + //xi = xinew; + //if ( ((xinew[0] + xinew[1] + xinew[2]) < 1) && + // (xinew[0] > 0) && (xinew[1] > 0) && (xinew[2] > 0) ) { + value2 = ma::getInsideness(mesh, apf::getMeshEntity(elems[i]), xinew); if ( value2 > bestValue) { bestValue = value2; elemNum = i; xi = xinew; + //std::cout<<" converged fval "<getCoordinateField(), me); + + for (int j = 0; j < iter; j++) { + apf::getJacobianInv(me, xin, Jinv); + //mapLocalToGlobal(me, xin, xyz); + apf::getVector(e, xin, xyz); + apf::Vector3 fval = (xyz-point); + + xinew = xin - transpose(Jinv)*fval; + + if ( fval.getLength() < tol ) { + xi = xinew; + //std::cout<<" fval "< elems(n); - for (int i = 0; i < n; i++) + apf::NewArray elemsF(n); + for (int i = 0; i < n; i++) { elems[i] = apf::createMeshElement(mesh, oldElements[i]); + for (int j = 0; j < mesh->countFields(); j++) { + apf::Field* fCrd = mesh->getField(j); + if (apf::getName(fCrd) == "CoordField") { + elemsF[i] = apf::createElement(fCrd, elems[i]); + } + } + } + // For a point outside the cavity getBestElement routine // returns an extrapolated xi // should perform a check for negative xi values @@ -273,11 +323,40 @@ namespace crv { &elems[0], xyz, xi); PCU_ALWAYS_ASSERT(entNum >= 0); + apf::Vector3 xiBetter = xi; + bool betterXi = getBetterXiApproximate(mesh, elems[entNum],xyz, xiBetter); + if (betterXi) xi = xiBetter; +/* + apf::MeshElement* meshElemP = apf::createMeshElement( + mesh, oldElements[entNum]); + apf::Vector3 xxyy, xxzz, val2; + apf::mapLocalToGlobal(meshElemP, xi, xxyy); + const char* namef = apf::getName(f); + std::cout<<" xi "<getModelType(mesh->toModel(newEnt))< val1; + val1.allocate(apf::countComponents(f)); + if (apf::getValueType(f) == apf::VECTOR) { + apf::mapLocalToGlobal(meshElemP, xi, xxzz); + apf::getVector(elemPC, xi, val2); + + apf::getComponents(elemPC, xi, &(val1[0])); + std::cout<<" Allocated values ("<< val2[0]-xxzz[0]<<","<< + val2[1]-xxzz[1]<<","< val; val.allocate(apf::countComponents(field)); apf::getComponents(elemP, xi, &(val[0])); + apf::setComponents(f, newEnt, nodeNum, &(val[0])); for (int i = 0; i < n; i++) @@ -300,27 +379,43 @@ namespace crv { int td = apf::Mesh::typeDimension[type]; if (td != d) continue; - int order = shape->getOrder(); int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); if (!ne) continue; apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); + apf::Element* e = apf::createElement(mesh->getCoordinateField(), me); apf::Vector3 point, xinew; - int n = shape->getEntityShape(type)->countNodes(); for (int j = 0; j < ne; j++) { shape->getNodeXi(type, j, xinew); - apf::mapLocalToGlobal(me, xinew, point); + apf::getVector(e, xinew, point); + //apf::mapLocalToGlobal(me, xinew, point); setInterpolatingFieldValues(mesh, oldElements, newEntities[i], j, point, f); } + apf::destroyMeshElement(me); + } + } + // convert interpolating values to control points in the + // decreasing order of entity dimension + for (int d = 3; d <= 1; d++) { + for(size_t i = 0; i < newEntities.getSize(); i++) { + int type = mesh->getType(newEntities[i]); + int td = apf::Mesh::typeDimension[type]; + if (td != d) continue; + + int order = shape->getOrder(); + int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); + + if (!ne) continue; + + int n = shape->getEntityShape(type)->countNodes(); apf::NewArray c; crv::getBezierTransformationCoefficients(order, apf::Mesh::simplexTypes[td], c); crv::convertInterpolationFieldPoints(newEntities[i], f, n, ne, c); - apf::destroyMeshElement(me); } } } diff --git a/crv/crvCurveMesh.cc b/crv/crvCurveMesh.cc index 966706b30..d63aec45e 100644 --- a/crv/crvCurveMesh.cc +++ b/crv/crvCurveMesh.cc @@ -193,12 +193,12 @@ bool BezierCurver::run() } convertInterpolatingToBezier(); - +/* if( m_mesh->getDimension() >= 2 && m_order == 2){ ma::Input* shapeFixer = configureShapeCorrection(m_mesh); crv::adapt(shapeFixer); } - +*/ m_mesh->acceptChanges(); m_mesh->verify(); return true; diff --git a/crv/crvShape.cc b/crv/crvShape.cc index 33f421292..1658ea7c8 100644 --- a/crv/crvShape.cc +++ b/crv/crvShape.cc @@ -751,6 +751,7 @@ static void swapInvalidEdges(Adapt* a) static void repositionInvalidEdges(Adapt* a) { + return; double t0 = PCU_Time(); EdgeReshaper es(a); ma::applyOperator(a,&es); @@ -836,17 +837,17 @@ void fixCrvElementShapes(Adapt* a) if ( ! count) break; prev_count = count; - fixLargeAngles(a); // update this - /* int numOpSuccess = fixLargeAngles(a); // update this */ - /* PCU_Add_Ints(&numOpSuccess,1); */ - /* if (PCU_Comm_Self() == 0) */ - /* lion_oprint(1,"==> %d large angle fix operations succeeded.\n", numOpSuccess); */ - markCrvBadQuality(a); + //fixLargeAngles(a); // update this + //int numOpSuccess = fixLargeAngles(a); // update this */ + //PCU_Add_Ints(&numOpSuccess,1); + //if (PCU_Comm_Self() == 0) + // lion_oprint(1,"==> %d large angle fix operations succeeded.\n", numOpSuccess); + //markCrvBadQuality(a); fixShortEdgeElements(a); // update this - /* int numEdgeRemoved = fixShortEdgeElements(a); // update this */ - /* PCU_Add_Ints(&numEdgeRemoved,1); */ - /* if (PCU_Comm_Self() == 0) */ - /* lion_oprint(1,"==> %d edges removal operations succeeded.\n", numEdgeRemoved); */ + //int numEdgeRemoved = fixShortEdgeElements(a); // update this */ + //PCU_Add_Ints(&numEdgeRemoved,1); + //if (PCU_Comm_Self() == 0) + // lion_oprint(1,"==> %d edges removal operations succeeded.\n", numEdgeRemoved); count = markCrvBadQuality(a); ++i; } while(count < prev_count && i < 6); // the second conditions is to make sure this does not take long diff --git a/ma/maShape.cc b/ma/maShape.cc index 458c9f4eb..dbf42b3aa 100644 --- a/ma/maShape.cc +++ b/ma/maShape.cc @@ -814,6 +814,7 @@ void fixElementShapes(Adapt* a) return; double t0 = PCU_Time(); int count = markBadQuality(a); + printf(" number of bad shape inside fixElementShape %d \n", count); int originalCount = count; int prev_count; double time; diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index bdfcbcea9..28a4aef98 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -229,7 +229,7 @@ void CavityTransfer::onCavity( double val = apf::getScalar(field, newEntities[i], 0); //double val = apf::measure(mesh, newEntities[i]); //apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); - if (val == 0) continue; + if (sumNew == 0) continue; apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); } } From 251213b74280f590f165308c10c436414ee3c8f1 Mon Sep 17 00:00:00 2001 From: Avinash Date: Mon, 3 Aug 2020 07:12:11 -0400 Subject: [PATCH 28/36] debugging field values after solutionTransfer + commented second order logSizeField --- crv/crvAdapt.cc | 40 +- crv/crvBezierSolutionTransfer.cc | 780 +++++++++++++++++-------------- crv/crvShape.cc | 41 +- crv/crvShapeHandler.cc | 8 +- ma/maSize.cc | 24 + 5 files changed, 547 insertions(+), 346 deletions(-) diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index 330952fd6..53c009777 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -17,8 +17,46 @@ #include #include #include - +#include +#include namespace crv { +//DEBUGGING +static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) +{ + apf::MeshEntity* e; + apf::MeshIterator* it = mesh->begin(3); + apf::Vector3 val1, val2, diff; + bool isSame = true; + + while ( (e = mesh->iterate(it)) ) { + apf::MeshElement* me = apf::createMeshElement(mesh, e); + apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); + apf::Element* e2; + for (int j = 0; j < mesh->countFields(); j++) { + apf::Field* fCrd = mesh->getField(j); + const char* namef = apf::getName(fCrd); + if (strcmp(namef,"CoordField") == 0) + e2 = apf::createElement(fCrd, me); + } + apf::getVector(e1, xi, val1); + apf::getVector(e2, xi, val2); + + diff = val1-val2; + //std::cout<<" values "<end(it); + + std::cout<<"---------------------------------------"< namespace crv { - static void getVertParams(apf::Mesh* mesh, int ptype, - apf::MeshEntity** parentVerts, - apf::NewArray& midEdgeVerts, - apf::MeshEntity* e, - apf::Vector3 params[4]) - { - int npv = apf::Mesh::adjacentCount[ptype][0]; - int ne = apf::Mesh::adjacentCount[ptype][1]; - apf::Downward verts; - - int nv = mesh->getDownward(e,0,verts); - // first check verts - for (int v = 0; v < nv; ++v){ - bool vert = false; - for (int i = 0; i < npv; ++i){ - if(verts[v] == parentVerts[i]){ - params[v] = elem_vert_xi[ptype][i]; - vert = true; - break; - } +static bool compareFields(apf::Mesh* mesh, + apf::MeshEntity* e, apf::Vector3 xi) +{ + apf::Vector3 val1, val2, diff; + bool isSame = true; + + apf::MeshElement* me = apf::createMeshElement(mesh, e); + apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); + apf::Element* e2; + for (int j = 0; j < mesh->countFields(); j++) { + apf::Field* fCrd = mesh->getField(j); + const char* namef = apf::getName(fCrd); + if (strcmp(namef,"CoordField") == 0) + e2 = apf::createElement(fCrd, me); + } + apf::getVector(e1, xi, val1); + apf::getVector(e2, xi, val2); + + diff = val1-val2; + std::cout<<" values "<countNodesOn(apf::Mesh::EDGE); + m->getDownward(edge,0,verts); + /* + apf::MeshElement* me1 = apf::createMeshElement(m, verts[0]); + apf::MeshElement* me2 = apf::createMeshElement(m, verts[1]); + apf::Element* e1 = apf::createElement(f, me1); + apf::Element* e2 = apf::createElement(f, me2); + */ + apf::NewArray value1, value2, value; + value.allocate(apf::countComponents(f)); + value1.allocate(apf::countComponents(f)); + value2.allocate(apf::countComponents(f)); + apf::getComponents(f, verts[0], 0, &(value1[0])); + apf::getComponents(f, verts[1], 0, &(value2[0])); + + m->getPoint(verts[0],0,points[0]); + m->getPoint(verts[1],0,points[1]); + for (int j = 0; j < ni; ++j){ + double t = (1.+j)/(1.+ni); + for (int k = 0; k < apf::countComponents(f); k++) + value[k] = value1[k]*(1.-t) + value2[k]*t; + apf::setComponents(f, edge, j, &(value[0])); + } + //apf::destroyElement(e); + //apf::destroyMeshElement(me); +} + +static void repositionInteriorFieldWithBlended(ma::Mesh* m, + apf::Field* f, ma::Entity* e) +{ + apf::FieldShape * fs = apf::getShape(f); + int order = fs->getOrder(); + int typeDim = apf::Mesh::typeDimension[m->getType(e)]; + int type = apf::Mesh::simplexTypes[typeDim]; + + if(!fs->hasNodesIn(typeDim) || getBlendingOrder(type)) + return; + + int n = fs->getEntityShape(type)->countNodes(); + int ne = fs->countNodesOn(type); + apf::NewArray c; + crv::getInternalBezierTransformationCoefficients(m,order,1,type,c); + crv::convertInterpolationFieldPoints(e,f,n-ne,ne,c); + +} + + +static void getVertParams(apf::Mesh* mesh, int ptype, + apf::MeshEntity** parentVerts, + apf::NewArray& midEdgeVerts, + apf::MeshEntity* e, + apf::Vector3 params[4]) +{ + int npv = apf::Mesh::adjacentCount[ptype][0]; + int ne = apf::Mesh::adjacentCount[ptype][1]; + apf::Downward verts; + + int nv = mesh->getDownward(e,0,verts); + // first check verts + for (int v = 0; v < nv; ++v){ + bool vert = false; + for (int i = 0; i < npv; ++i){ + if(verts[v] == parentVerts[i]){ + params[v] = elem_vert_xi[ptype][i]; + vert = true; + break; } + } - if(!vert){ - for (int i = 0; i < ne; ++i){ - if( verts[v] == midEdgeVerts[i] ){ - params[v] = elem_edge_xi[ptype][i]; - break; - } - } + if(!vert){ + for (int i = 0; i < ne; ++i){ + if( verts[v] == midEdgeVerts[i] ){ + params[v] = elem_edge_xi[ptype][i]; + break; + } } } } +} - static int getBestElement(int n, - apf::Mesh* mesh, - apf::MeshElement **elems, - apf::Vector3 point, apf::Vector3 &xi) - { - int iter = 1000; - double tol = 1e-4; - - int elemNum = 0; - double value2, bestValue = -DBL_MAX; - - apf::Vector3 initialGuess = apf::Vector3(1./4., 1./4., 1./4.); - apf::Vector3 xinew, xyz; - apf::Matrix3x3 Jinv; - apf::Vector3 fval; - - //printf(" inside find best element loop %d \n", n); - for (int i = 0; i < n; i++) { - apf::Vector3 xin = initialGuess; - apf::Element* e = apf::createElement(mesh->getCoordinateField(), elems[i]); - - for (int j = 0; j < iter; j++) { - apf::getJacobianInv(elems[i], xin, Jinv); - //mapLocalToGlobal(elems[i], xin, xyz); - apf::getVector(e, xin, xyz); - fval = (xyz-point); - - xinew = xin - transpose(Jinv)*fval; - - if ( (xinew-xin).getLength() < tol ) { - //elemNum = i; - //xi = xinew; - //if ( ((xinew[0] + xinew[1] + xinew[2]) < 1) && - // (xinew[0] > 0) && (xinew[1] > 0) && (xinew[2] > 0) ) { - - value2 = ma::getInsideness(mesh, apf::getMeshEntity(elems[i]), xinew); - if ( value2 > bestValue) { - bestValue = value2; - elemNum = i; - xi = xinew; - //std::cout<<" converged fval "<getCoordinateField(), me); + //printf(" inside find best element loop %d \n", n); + for (int i = 0; i < n; i++) { + apf::Vector3 xin = initialGuess; + apf::Element* e = apf::createElement(mesh->getCoordinateField(), elems[i]); for (int j = 0; j < iter; j++) { - apf::getJacobianInv(me, xin, Jinv); - //mapLocalToGlobal(me, xin, xyz); + apf::getJacobianInv(elems[i], xin, Jinv); + //mapLocalToGlobal(elems[i], xin, xyz); apf::getVector(e, xin, xyz); - apf::Vector3 fval = (xyz-point); + fval = (xyz-point); xinew = xin - transpose(Jinv)*fval; - if ( fval.getLength() < tol ) { - xi = xinew; - //std::cout<<" fval "< bestValue) { + bestValue = value2; + elemNum = i; + xi = xinew; + //std::cout<<" converged fval "<getCoordinateField(), me); + + for (int j = 0; j < iter; j++) { + apf::getJacobianInv(me, xin, Jinv); + //mapLocalToGlobal(me, xin, xyz); + apf::getVector(e, xin, xyz); + apf::Vector3 fval = (xyz-point); + + xinew = xin - transpose(Jinv)*fval; + + if ( fval.getLength() < tol ) { + xi = xinew; + //std::cout<<" fval "<refine; + mesh = adapt->mesh; + shape = apf::getShape(f); + int P = mesh->getShape()->getOrder(); + for (int d = 1; d <= mesh->getDimension(); d++) { + int type = apf::Mesh::simplexTypes[d]; + if (!mesh->getShape()->hasNodesIn(d)) + continue; + int n = mesh->getShape()->getEntityShape(type)->countNodes(); + mth::Matrix A(n,n); + Ai[d].resize(n,n); + getBezierTransformationMatrix(type, P, + A, elem_vert_xi[type]); + invertMatrixWithPLU(getNumControlPoints(type,P), A, Ai[d]); + } + } + virtual bool hasNodesOn(int dimension) { - f = field; - refine = adapt->refine; - mesh = adapt->mesh; - shape = apf::getShape(f); - int P = mesh->getShape()->getOrder(); - for (int d = 1; d <= mesh->getDimension(); d++) { - int type = apf::Mesh::simplexTypes[d]; - if (!mesh->getShape()->hasNodesIn(d)) - continue; - int n = mesh->getShape()->getEntityShape(type)->countNodes(); - mth::Matrix A(n,n); - Ai[d].resize(n,n); - getBezierTransformationMatrix(type, P, - A, elem_vert_xi[type]); - invertMatrixWithPLU(getNumControlPoints(type,P), A, Ai[d]); - } + return others.hasNodesOn(dimension); } - virtual bool hasNodesOn(int dimension) - { - return others.hasNodesOn(dimension); - } - virtual const char* getTransferFieldName() { - return apf::getName(f); - } + virtual const char* getTransferFieldName() { + return apf::getName(f); + } - virtual void onVertex( - apf::MeshElement* parent, - ma::Vector const& xi, - ma::Entity* vert) - { - apf::Element* e = apf::createElement(f, parent); - apf::NewArray value; - value.allocate(apf::countComponents(f)); - apf::getComponents(e, xi, &(value[0])); - apf::setComponents(f, vert, 0, &(value[0])); - apf::destroyElement(e); - } + virtual void onVertex( + apf::MeshElement* parent, + ma::Vector const& xi, + ma::Entity* vert) + { + apf::Element* e = apf::createElement(f, parent); + apf::NewArray value; + value.allocate(apf::countComponents(f)); + apf::getComponents(e, xi, &(value[0])); + apf::setComponents(f, vert, 0, &(value[0])); + apf::destroyElement(e); + } - virtual void onRefine( - ma::Entity* parent, - ma::EntityArray& newEntities) - { - int P = shape->getOrder(); - int parentType = mesh->getType(parent); - apf::Downward parentVerts, parentEdges; + virtual void onRefine( + ma::Entity* parent, + ma::EntityArray& newEntities) + { + int P = shape->getOrder(); + int parentType = mesh->getType(parent); + apf::Downward parentVerts, parentEdges; - mesh->getDownward(parent, 0, parentVerts); - mesh->getDownward(parent, 1, parentEdges); + mesh->getDownward(parent, 0, parentVerts); + mesh->getDownward(parent, 1, parentEdges); - int ne = apf::Mesh::adjacentCount[parentType][1]; + int ne = apf::Mesh::adjacentCount[parentType][1]; - //std::cout<<" field name "< midEdgeVerts(ne); + apf::NewArray midEdgeVerts(ne); + bool useLinear = false; - for (int i = 0; i < ne; ++i){ - if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) - midEdgeVerts[i] = ma::findSplitVert(refine,parentEdges[i]); - else - midEdgeVerts[i] = 0; + for (int i = 0; i < ne; ++i){ + if ( ma::getFlag(adapt,parentEdges[i],ma::SPLIT) ) + midEdgeVerts[i] = ma::findSplitVert(refine,parentEdges[i]); + else + midEdgeVerts[i] = 0; + if ( ma::getFlag(adapt,parentEdges[i], ma::BAD_QUALITY) ) { + useLinear = true; } + } - int np = shape->getEntityShape(parentType)->countNodes(); + int np = shape->getEntityShape(parentType)->countNodes(); - apf::Element* elem = apf::createElement(f, parent); - apf::NewArray Vnodes; - apf::NewArray Snodes; + apf::Element* elem = apf::createElement(f, parent); + apf::NewArray Vnodes; + apf::NewArray Snodes; - if (apf::getValueType(f) == apf::VECTOR) { - apf::getVectorNodes(elem,Vnodes); - } - else if (apf::getValueType(f) == apf::SCALAR) { - apf::getScalarNodes(elem,Snodes); + apf::NewArray Vcoord; + apf::Element* elemC = apf::createElement(mesh->getCoordinateField(), parent); + apf::getVectorNodes(elem,Vcoord); + + + if (apf::getValueType(f) == apf::VECTOR) { + apf::getVectorNodes(elem,Vnodes); + for (int kk = 0; kk < mesh->countFields(); kk++) { + apf::Field* crdf = mesh->getField(kk); + const char* namefc = apf::getName(crdf); + if (strcmp(namefc, "CoordField") == 0) { + for(int i = 0; i < np; i++) { + //std::cout<<" node diff "<hasNodesIn(d)) continue; - for (size_t i = 0; i < newEntities.getSize(); i++) { - int childType = mesh->getType(newEntities[i]); - if (d != apf::Mesh::typeDimension[childType]) - continue; + for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++) { + if (!shape->hasNodesIn(d)) continue; + for (size_t i = 0; i < newEntities.getSize(); i++) { + int childType = mesh->getType(newEntities[i]); + if (d != apf::Mesh::typeDimension[childType]) + continue; + int n = shape->getEntityShape(childType)->countNodes(); + int ni = shape->countNodesOn(childType); + + if (useLinear && !isBoundaryEntity(mesh, newEntities[i])) { + if (childType == apf::Mesh::EDGE) { + setLinearEdgeFieldPoints(mesh, f, newEntities[i]); + } else { + for (int j = 0; j < ni; j++) { + apf::NewArray val; + val.allocate(apf::countComponents(f)); + for (int k = 0; k < apf::countComponents(f); k++) + val[k] = 0.; + apf::setComponents(f, newEntities[i], j, &(val[0])); + } + repositionInteriorFieldWithBlended(mesh,f,newEntities[i]); + } + } + else { //int n = getNumControlPoints(childType,P); - int n = shape->getEntityShape(childType)->countNodes(); - int ni = shape->countNodesOn(childType); apf::Vector3 vp[4]; getVertParams(mesh, parentType,parentVerts, midEdgeVerts,newEntities[i], vp); @@ -255,9 +373,27 @@ namespace crv { if (apf::getValueType(f) == apf::VECTOR) { apf::Vector3 Vvalue(0,0,0); - for (int k = 0; k < np; ++k) + apf::Vector3 vcoordV(0,0,0); + for (int k = 0; k < np; ++k) { Vvalue += Vnodes[k]*B(j+n-ni,k); + } + //std::cout<<" field ("<getCoordinateField(), + newEntities[i], j, vcoordV); + for (int kk = 0; kk < mesh->countFields(); kk++) { + apf::Field* crdf = mesh->getField(kk); + const char* namefc = apf::getName(crdf); + if (strcmp(namefc, "CoordField") == 0) { + for (int kkk = 0; kkk < np; kkk++) { + //std::cout<<" inside allocation "<< + // Vvalue-vcoordV<= 1; --d){ - for (int d = 1; d <= apf::Mesh::typeDimension[parentType]; d++){ - if (!shape->hasNodesIn(d)) continue; - - for(size_t i = 0; i getType(newEntities[i]); - int td = apf::Mesh::typeDimension[type]; - if (td != d) continue; - - int ni = shape->countNodesOn(apf::Mesh::simplexTypes[td]); - - if (!ni) continue; - int n = shape->getEntityShape(type)->countNodes(); - - apf::NewArray c; - crv::getBezierTransformationCoefficients(P, - apf::Mesh::simplexTypes[td], c); - crv::convertInterpolationFieldPoints(newEntities[i], - f, n, ni, c); - } - } - */ - apf::destroyElement(elem); } - - void setInterpolatingFieldValues(apf::Mesh* mesh, - ma::EntityArray &oldElements, - apf::MeshEntity* newEnt, int nodeNum, - apf::Vector3 xyz, apf::Field* field) - { - apf::Vector3 xi; - int entNum = -1; - int n = oldElements.getSize(); - - apf::NewArray elems(n); - apf::NewArray elemsF(n); - for (int i = 0; i < n; i++) { - elems[i] = apf::createMeshElement(mesh, oldElements[i]); - for (int j = 0; j < mesh->countFields(); j++) { - apf::Field* fCrd = mesh->getField(j); - if (apf::getName(fCrd) == "CoordField") { - elemsF[i] = apf::createElement(fCrd, elems[i]); - } - } - } - - // For a point outside the cavity getBestElement routine - // returns an extrapolated xi - // should perform a check for negative xi values - entNum = getBestElement(n, mesh, - &elems[0], xyz, xi); - - PCU_ALWAYS_ASSERT(entNum >= 0); - apf::Vector3 xiBetter = xi; - bool betterXi = getBetterXiApproximate(mesh, elems[entNum],xyz, xiBetter); - if (betterXi) xi = xiBetter; -/* - apf::MeshElement* meshElemP = apf::createMeshElement( - mesh, oldElements[entNum]); - apf::Vector3 xxyy, xxzz, val2; - apf::mapLocalToGlobal(meshElemP, xi, xxyy); - const char* namef = apf::getName(f); - std::cout<<" xi "<getModelType(mesh->toModel(newEnt))< val1; - val1.allocate(apf::countComponents(f)); - if (apf::getValueType(f) == apf::VECTOR) { - apf::mapLocalToGlobal(meshElemP, xi, xxzz); - apf::getVector(elemPC, xi, val2); - - apf::getComponents(elemPC, xi, &(val1[0])); - std::cout<<" Allocated values ("<< val2[0]-xxzz[0]<<","<< - val2[1]-xxzz[1]<<","< val; - val.allocate(apf::countComponents(field)); - apf::getComponents(elemP, xi, &(val[0])); - - apf::setComponents(f, newEnt, nodeNum, &(val[0])); - - for (int i = 0; i < n; i++) - apf::destroyMeshElement(elems[i]); - - apf::destroyElement(elemP); + } + void setInterpolatingFieldValues(apf::Mesh* mesh, + ma::EntityArray &oldElements, + apf::MeshEntity* newEnt, int nodeNum, + apf::Vector3 xyz, apf::Field* field1) + { + apf::Vector3 xi; + int entNum = -1; + int n = oldElements.getSize(); + apf::FieldShape* fshape = apf::getShape(field1); + + apf::NewArray elems(n); + apf::NewArray elemsF(n); + for (int i = 0; i < n; i++) { + elems[i] = apf::createMeshElement(mesh, oldElements[i]); } - virtual void onCavity( - ma::EntityArray& oldElements, - ma::EntityArray& newEntities) - { - - if (getDimension(mesh, oldElements[0]) < 3) - return; - - for (int d = 1; d < 4; d++) { - for(size_t i = 0; i getType(newEntities[i]); - int td = apf::Mesh::typeDimension[type]; - if (td != d) continue; - - int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); + // For a point outside the cavity getBestElement routine + // returns an extrapolated xi + // should perform a check for negative xi values + entNum = getBestElement(n, mesh, + &elems[0], xyz, xi); + + PCU_ALWAYS_ASSERT(entNum >= 0); + apf::Vector3 xiBetter = xi; + //bool betterXi = getBetterXiApproximate(mesh, elems[entNum],xyz, xiBetter); + //if (betterXi) xi = xiBetter; + + apf::MeshElement* meshElemP = apf::createMeshElement( + mesh, oldElements[entNum]); + apf::Vector3 xxyy, xxzz, val2; + apf::mapLocalToGlobal(meshElemP, xi, xxyy); + const char* namef = apf::getName(field1); + apf::NewArray shpVal; + apf::NewArray shpVal1; + apf::Element* elemFld = apf::createElement( + mesh->getCoordinateField(), meshElemP); + apf::getShapeValues(elemFld, xi, shpVal); + + int parentType = mesh->getType(oldElements[entNum]); + int np = fshape->getEntityShape(parentType)->countNodes(); + + apf::Element* elemP = apf::createElement(field1, oldElements[entNum]); + apf::NewArray val; + val.allocate(apf::countComponents(field1)); + apf::getComponents(elemP, xi, &(val[0])); + + apf::setComponents(field1, newEnt, nodeNum, &(val[0])); + for (int i = 0; i < n; i++) + apf::destroyMeshElement(elems[i]); + + apf::destroyElement(elemP); + } - if (!ne) continue; + virtual void onCavity( + ma::EntityArray& oldElements, + ma::EntityArray& newEntities) + { - apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); - apf::Element* e = apf::createElement(mesh->getCoordinateField(), me); - apf::Vector3 point, xinew; + if (getDimension(mesh, oldElements[0]) < 3) + return; - for (int j = 0; j < ne; j++) { - shape->getNodeXi(type, j, xinew); - apf::getVector(e, xinew, point); - //apf::mapLocalToGlobal(me, xinew, point); - setInterpolatingFieldValues(mesh, oldElements, - newEntities[i], j, point, f); - } - apf::destroyMeshElement(me); - } - } - // convert interpolating values to control points in the - // decreasing order of entity dimension - for (int d = 3; d <= 1; d++) { - for(size_t i = 0; i < newEntities.getSize(); i++) { - int type = mesh->getType(newEntities[i]); - int td = apf::Mesh::typeDimension[type]; - if (td != d) continue; + for (int d = 1; d < 4; d++) { + for(size_t i = 0; i getType(newEntities[i]); + int td = apf::Mesh::typeDimension[type]; + if (td != d) continue; - int order = shape->getOrder(); - int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); + int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); - if (!ne) continue; + if (!ne) continue; - int n = shape->getEntityShape(type)->countNodes(); - apf::NewArray c; - crv::getBezierTransformationCoefficients(order, - apf::Mesh::simplexTypes[td], c); + apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); + apf::Element* e = apf::createElement(mesh->getCoordinateField(), me); + apf::Vector3 point, xinew; - crv::convertInterpolationFieldPoints(newEntities[i], f, n, ne, c); + for (int j = 0; j < ne; j++) { + shape->getNodeXi(type, j, xinew); + apf::getVector(e, xinew, point); + //apf::mapLocalToGlobal(me, xinew, point); + setInterpolatingFieldValues(mesh, oldElements, + newEntities[i], j, point, f); } + + apf::destroyMeshElement(me); } } + // convert interpolating values to control points in the + // decreasing order of entity dimension + for (int d = 3; d >= 1; d--) { + for(size_t i = 0; i < newEntities.getSize(); i++) { + int type = mesh->getType(newEntities[i]); + int td = apf::Mesh::typeDimension[type]; + if (td != d) continue; - protected: - apf::Field* f; - ma::Mesh* mesh; - ma::Refine* refine; - apf::FieldShape* shape; - mth::Matrix Ai[4]; - }; - - static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* a) - { - return new CrvBezierSolutionTransfer(f, a); - } + int order = shape->getOrder(); + int ne = shape->countNodesOn(apf::Mesh::simplexTypes[td]); - ma::SolutionTransfer* setBezierSolutionTransfers( - const std::vector& fields, Adapt* a) - { - ma::SolutionTransfers *st = dynamic_cast(a->solutionTransfer); - if (!st) - st = new ma::SolutionTransfers(); + if (!ne) continue; - for (std::size_t i = 0; i < fields.size(); i++) { - st->add(createBezierSolutionTransfer(fields[i], a)); - } - std::vector trans = st->transfers; - for (std::size_t i = 0; i < trans.size(); i++) { - const char* name = trans[i]->getTransferFieldName(); - std::cout<<" field name added "<< name<getEntityShape(type)->countNodes(); + apf::NewArray c; + crv::getBezierTransformationCoefficients(order, + apf::Mesh::simplexTypes[td], c); + + crv::convertInterpolationFieldPoints(newEntities[i], f, n, ne, c); } - return st; } + } + protected: + apf::Field* f; + ma::Mesh* mesh; + ma::Refine* refine; + apf::FieldShape* shape; + mth::Matrix Ai[4]; +}; + +static ma::SolutionTransfer* createBezierSolutionTransfer(apf::Field* f, Adapt* a) +{ + return new CrvBezierSolutionTransfer(f, a); +} + +ma::SolutionTransfer* setBezierSolutionTransfers( + const std::vector& fields, Adapt* a) +{ + ma::SolutionTransfers *st = dynamic_cast(a->solutionTransfer); + if (!st) + st = new ma::SolutionTransfers(); + + for (std::size_t i = 0; i < fields.size(); i++) { + st->add(createBezierSolutionTransfer(fields[i], a)); } + std::vector trans = st->transfers; + for (std::size_t i = 0; i < trans.size(); i++) { + const char* name = trans[i]->getTransferFieldName(); + std::cout<<" field name added "<< name< #include #include - +#include /* This is similar to maShape.cc, conceptually, but different enough * that some duplicate code makes sense */ namespace crv { +static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) +{ + apf::MeshEntity* e; + apf::MeshIterator* it = mesh->begin(3); + apf::Vector3 val1, val2, diff; + bool isSame = true; + + while ( (e = mesh->iterate(it)) ) { + apf::MeshElement* me = apf::createMeshElement(mesh, e); + apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); + apf::Element* e2; + for (int j = 0; j < mesh->countFields(); j++) { + apf::Field* fCrd = mesh->getField(j); + const char* namef = apf::getName(fCrd); + if (strcmp(namef,"CoordField") == 0) + e2 = apf::createElement(fCrd, me); + } + apf::getVector(e1, xi, val1); + apf::getVector(e2, xi, val2); + + diff = val1-val2; + //std::cout<<" diff "<end(it); + //std::cout<<" ================================" <getModelType(m->toModel(e)) < m->getDimension(); diff --git a/crv/crvShapeHandler.cc b/crv/crvShapeHandler.cc index 76e37c435..7d4bf41ca 100644 --- a/crv/crvShapeHandler.cc +++ b/crv/crvShapeHandler.cc @@ -19,7 +19,7 @@ #include #include #include - +#include namespace crv { static double measureLinearTriArea(ma::Mesh* m, ma::Entity* tri) @@ -189,6 +189,7 @@ class BezierTransfer : public ma::SolutionTransfer int n = getNumControlPoints(childType,P); apf::Vector3 vp[4]; getVertParams(parentType,parentVerts,midEdgeVerts,newEntities[i],vp); + if(getBlendingOrder(childType) > 0){ apf::NewArray childXi(n); collectNodeXi(parentType,childType,P,vp,childXi); @@ -203,8 +204,11 @@ class BezierTransfer : public ma::SolutionTransfer mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); for (int j = 0; j < ni; ++j){ apf::Vector3 point(0,0,0); - for (int k = 0; k < np; ++k) + for (int k = 0; k < np; ++k) { point += nodes[k]*B(j+n-ni,k); + } + //std::cout<<" shape ("<setPoint(newEntities[i],j,point); } } diff --git a/ma/maSize.cc b/ma/maSize.cc index c0fab50c7..64b64d199 100644 --- a/ma/maSize.cc +++ b/ma/maSize.cc @@ -429,8 +429,11 @@ struct LogAnisoSizeField : public MetricSizeField logMEval(f) { mesh = m; + // ADDED logMField = apf::createUserField(m, "ma_logM", apf::MATRIX, apf::getLagrange(1), &logMEval); + //logMField = apf::createUserField(m, "ma_logM", apf::MATRIX, + // apf::getLagrange(2), &logMEval); } ~LogAnisoSizeField() { @@ -439,7 +442,9 @@ struct LogAnisoSizeField : public MetricSizeField void init(Mesh* m, apf::Field* sizes, apf::Field* frames) { mesh = m; + //ADDED logMField = apf::createField(m, "ma_logM", apf::MATRIX, apf::getLagrange(1)); + //logMField = apf::createField(m, "ma_logM", apf::MATRIX, apf::getLagrange(2)); Entity* v; Iterator* it = m->begin(0); while ( (v = m->iterate(it)) ) { @@ -453,6 +458,25 @@ struct LogAnisoSizeField : public MetricSizeField 0 , 0 , s[2]); apf::setMatrix(logMField, v, 0, f * S * transpose(f)); } + m->end(it); + + //ADDED + /* + it = m->begin(1); + while ( (v = m->iterate(it)) ) { + Vector h; + Matrix f; + apf::getVector(sizes, v, 0, h); + apf::getMatrix(frames, v, 0, f); + Vector s(log(1/h[0]/h[0]), log(1/h[1]/h[1]), log(1/h[2]/h[2])); + Matrix S(s[0], 0 , 0, + 0 , s[1], 0, + 0 , 0 , s[2]); + apf::setMatrix(logMField, v, 0, f * S * transpose(f)); + } + m->end(it); + */ + } void getTransform( apf::MeshElement* me, From 1db872437bf8db964786de60d0bf872c2428a8d1 Mon Sep 17 00:00:00 2001 From: Avinash Date: Tue, 18 Aug 2020 12:53:46 -0400 Subject: [PATCH 29/36] clean code. --- crv/crvBezierSolutionTransfer.cc | 62 ++------------------------------ crv/crvShape.cc | 12 +++---- 2 files changed, 7 insertions(+), 67 deletions(-) diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 1453214d6..169466d97 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -27,6 +27,8 @@ #include namespace crv { + +// FOR DEBUGGING static bool compareFields(apf::Mesh* mesh, apf::MeshEntity* e, apf::Vector3 xi) { @@ -70,12 +72,6 @@ static void setLinearEdgeFieldPoints(ma::Mesh* m, apf::Field* f, int ni = fs->countNodesOn(apf::Mesh::EDGE); m->getDownward(edge,0,verts); - /* - apf::MeshElement* me1 = apf::createMeshElement(m, verts[0]); - apf::MeshElement* me2 = apf::createMeshElement(m, verts[1]); - apf::Element* e1 = apf::createElement(f, me1); - apf::Element* e2 = apf::createElement(f, me2); - */ apf::NewArray value1, value2, value; value.allocate(apf::countComponents(f)); value1.allocate(apf::countComponents(f)); @@ -91,8 +87,6 @@ static void setLinearEdgeFieldPoints(ma::Mesh* m, apf::Field* f, value[k] = value1[k]*(1.-t) + value2[k]*t; apf::setComponents(f, edge, j, &(value[0])); } - //apf::destroyElement(e); - //apf::destroyMeshElement(me); } static void repositionInteriorFieldWithBlended(ma::Mesh* m, @@ -164,14 +158,12 @@ static int getBestElement(int n, apf::Matrix3x3 Jinv; apf::Vector3 fval; - //printf(" inside find best element loop %d \n", n); for (int i = 0; i < n; i++) { apf::Vector3 xin = initialGuess; apf::Element* e = apf::createElement(mesh->getCoordinateField(), elems[i]); for (int j = 0; j < iter; j++) { apf::getJacobianInv(elems[i], xin, Jinv); - //mapLocalToGlobal(elems[i], xin, xyz); apf::getVector(e, xin, xyz); fval = (xyz-point); @@ -183,7 +175,6 @@ static int getBestElement(int n, bestValue = value2; elemNum = i; xi = xinew; - //std::cout<<" converged fval "<refine; @@ -292,8 +280,6 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer int ne = apf::Mesh::adjacentCount[parentType][1]; - //std::cout<<" field name "< midEdgeVerts(ne); bool useLinear = false; @@ -313,28 +299,13 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer apf::NewArray Vnodes; apf::NewArray Snodes; - apf::NewArray Vcoord; - apf::Element* elemC = apf::createElement(mesh->getCoordinateField(), parent); - apf::getVectorNodes(elem,Vcoord); - - if (apf::getValueType(f) == apf::VECTOR) { apf::getVectorNodes(elem,Vnodes); - for (int kk = 0; kk < mesh->countFields(); kk++) { - apf::Field* crdf = mesh->getField(kk); - const char* namefc = apf::getName(crdf); - if (strcmp(namefc, "CoordField") == 0) { - for(int i = 0; i < np; i++) { - //std::cout<<" node diff "<hasNodesIn(d)) continue; for (size_t i = 0; i < newEntities.getSize(); i++) { @@ -360,7 +331,6 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer } else { - //int n = getNumControlPoints(childType,P); apf::Vector3 vp[4]; getVertParams(mesh, parentType,parentVerts, midEdgeVerts,newEntities[i], vp); @@ -373,27 +343,11 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer if (apf::getValueType(f) == apf::VECTOR) { apf::Vector3 Vvalue(0,0,0); - apf::Vector3 vcoordV(0,0,0); for (int k = 0; k < np; ++k) { Vvalue += Vnodes[k]*B(j+n-ni,k); } - //std::cout<<" field ("<getCoordinateField(), - newEntities[i], j, vcoordV); - for (int kk = 0; kk < mesh->countFields(); kk++) { - apf::Field* crdf = mesh->getField(kk); - const char* namefc = apf::getName(crdf); - if (strcmp(namefc, "CoordField") == 0) { - for (int kkk = 0; kkk < np; kkk++) { - //std::cout<<" inside allocation "<< - // Vvalue-vcoordV<= 0); - apf::Vector3 xiBetter = xi; - //bool betterXi = getBetterXiApproximate(mesh, elems[entNum],xyz, xiBetter); - //if (betterXi) xi = xiBetter; apf::MeshElement* meshElemP = apf::createMeshElement( mesh, oldElements[entNum]); apf::Vector3 xxyy, xxzz, val2; apf::mapLocalToGlobal(meshElemP, xi, xxyy); - const char* namef = apf::getName(field1); - apf::NewArray shpVal; - apf::NewArray shpVal1; - apf::Element* elemFld = apf::createElement( - mesh->getCoordinateField(), meshElemP); - apf::getShapeValues(elemFld, xi, shpVal); int parentType = mesh->getType(oldElements[entNum]); int np = fshape->getEntityShape(parentType)->countNodes(); diff --git a/crv/crvShape.cc b/crv/crvShape.cc index 9cf432624..29e3666f0 100644 --- a/crv/crvShape.cc +++ b/crv/crvShape.cc @@ -25,23 +25,23 @@ namespace crv { -static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) +static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) { apf::MeshEntity* e; apf::MeshIterator* it = mesh->begin(3); apf::Vector3 val1, val2, diff; bool isSame = true; - while ( (e = mesh->iterate(it)) ) { - apf::MeshElement* me = apf::createMeshElement(mesh, e); + while ( (e = mesh->iterate(it)) ) { + apf::MeshElement* me = apf::createMeshElement(mesh, e); apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); - apf::Element* e2; + apf::Element* e2; for (int j = 0; j < mesh->countFields(); j++) { apf::Field* fCrd = mesh->getField(j); const char* namef = apf::getName(fCrd); if (strcmp(namef,"CoordField") == 0) e2 = apf::createElement(fCrd, me); - } + } apf::getVector(e1, xi, val1); apf::getVector(e2, xi, val2); @@ -57,8 +57,6 @@ static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) apf::destroyElement(e2); } mesh->end(it); - //std::cout<<" ================================" < Date: Fri, 28 Aug 2020 15:13:04 -0400 Subject: [PATCH 30/36] Code cleanup. --- crv/crvAdapt.cc | 37 ----------------- crv/crvBezierSolutionTransfer.cc | 71 +------------------------------- crv/crvShape.cc | 37 ----------------- 3 files changed, 2 insertions(+), 143 deletions(-) diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index 53c009777..86214aabc 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -20,43 +20,6 @@ #include #include namespace crv { -//DEBUGGING -static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) -{ - apf::MeshEntity* e; - apf::MeshIterator* it = mesh->begin(3); - apf::Vector3 val1, val2, diff; - bool isSame = true; - - while ( (e = mesh->iterate(it)) ) { - apf::MeshElement* me = apf::createMeshElement(mesh, e); - apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); - apf::Element* e2; - for (int j = 0; j < mesh->countFields(); j++) { - apf::Field* fCrd = mesh->getField(j); - const char* namef = apf::getName(fCrd); - if (strcmp(namef,"CoordField") == 0) - e2 = apf::createElement(fCrd, me); - } - apf::getVector(e1, xi, val1); - apf::getVector(e2, xi, val2); - - diff = val1-val2; - //std::cout<<" values "<end(it); - - std::cout<<"---------------------------------------"< namespace crv { - -// FOR DEBUGGING -static bool compareFields(apf::Mesh* mesh, - apf::MeshEntity* e, apf::Vector3 xi) -{ - apf::Vector3 val1, val2, diff; - bool isSame = true; - - apf::MeshElement* me = apf::createMeshElement(mesh, e); - apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); - apf::Element* e2; - for (int j = 0; j < mesh->countFields(); j++) { - apf::Field* fCrd = mesh->getField(j); - const char* namef = apf::getName(fCrd); - if (strcmp(namef,"CoordField") == 0) - e2 = apf::createElement(fCrd, me); - } - apf::getVector(e1, xi, val1); - apf::getVector(e2, xi, val2); - - diff = val1-val2; - std::cout<<" values "<getCoordinateField(), me); - - for (int j = 0; j < iter; j++) { - apf::getJacobianInv(me, xin, Jinv); - //mapLocalToGlobal(me, xin, xyz); - apf::getVector(e, xin, xyz); - apf::Vector3 fval = (xyz-point); - - xinew = xin - transpose(Jinv)*fval; - - if ( fval.getLength() < tol ) { - xi = xinew; - return true; - } - else - xin = xinew; - } - return false; -} - class CrvBezierSolutionTransfer : public ma::SolutionTransfer { public: @@ -369,7 +303,6 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer apf::Vector3 xi; int entNum = -1; int n = oldElements.getSize(); - apf::FieldShape* fshape = apf::getShape(field1); apf::NewArray elems(n); apf::NewArray elemsF(n); @@ -387,8 +320,8 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer apf::Vector3 xxyy, xxzz, val2; apf::mapLocalToGlobal(meshElemP, xi, xxyy); - int parentType = mesh->getType(oldElements[entNum]); - int np = fshape->getEntityShape(parentType)->countNodes(); + //int parentType = mesh->getType(oldElements[entNum]); + //int np = fshape->getEntityShape(parentType)->countNodes(); apf::Element* elemP = apf::createElement(field1, oldElements[entNum]); apf::NewArray val; diff --git a/crv/crvShape.cc b/crv/crvShape.cc index 29e3666f0..5c8d254b4 100644 --- a/crv/crvShape.cc +++ b/crv/crvShape.cc @@ -25,43 +25,6 @@ namespace crv { -static bool compareFields(apf::Mesh* mesh, apf::Vector3 xi) -{ - apf::MeshEntity* e; - apf::MeshIterator* it = mesh->begin(3); - apf::Vector3 val1, val2, diff; - bool isSame = true; - - while ( (e = mesh->iterate(it)) ) { - apf::MeshElement* me = apf::createMeshElement(mesh, e); - apf::Element* e1 = apf::createElement(mesh->getCoordinateField(), me); - apf::Element* e2; - for (int j = 0; j < mesh->countFields(); j++) { - apf::Field* fCrd = mesh->getField(j); - const char* namef = apf::getName(fCrd); - if (strcmp(namef,"CoordField") == 0) - e2 = apf::createElement(fCrd, me); - } - apf::getVector(e1, xi, val1); - apf::getVector(e2, xi, val2); - - diff = val1-val2; - //std::cout<<" diff "<end(it); - - return isSame; - -} - bool isBoundaryEntity(apf::Mesh* m, apf::MeshEntity* e) { return m->getModelType(m->toModel(e)) < m->getDimension(); From 721999b14c075b26e13b82504a7a833c30bb2d52 Mon Sep 17 00:00:00 2001 From: Avinash Date: Wed, 2 Sep 2020 01:26:20 -0400 Subject: [PATCH 31/36] cleanup --- crv/crvCurveMesh.cc | 5 ++--- crv/crvShape.cc | 23 +++++++++++------------ ma/CMakeLists.txt | 2 -- ma/maShape.cc | 1 - ma/maSize.cc | 24 ------------------------ 5 files changed, 13 insertions(+), 42 deletions(-) diff --git a/crv/crvCurveMesh.cc b/crv/crvCurveMesh.cc index d63aec45e..5cc39bbfb 100644 --- a/crv/crvCurveMesh.cc +++ b/crv/crvCurveMesh.cc @@ -150,7 +150,6 @@ void BezierCurver::convertInterpolatingToBezier() } // if we have a full representation, we need to place internal nodes on // triangles and tetrahedra - for(int d = 2; d <= md; ++d){ if(!fs->hasNodesIn(d) || getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; @@ -193,12 +192,12 @@ bool BezierCurver::run() } convertInterpolatingToBezier(); -/* + if( m_mesh->getDimension() >= 2 && m_order == 2){ ma::Input* shapeFixer = configureShapeCorrection(m_mesh); crv::adapt(shapeFixer); } -*/ + m_mesh->acceptChanges(); m_mesh->verify(); return true; diff --git a/crv/crvShape.cc b/crv/crvShape.cc index 5c8d254b4..33f421292 100644 --- a/crv/crvShape.cc +++ b/crv/crvShape.cc @@ -19,7 +19,7 @@ #include #include #include -#include + /* This is similar to maShape.cc, conceptually, but different enough * that some duplicate code makes sense */ @@ -751,7 +751,6 @@ static void swapInvalidEdges(Adapt* a) static void repositionInvalidEdges(Adapt* a) { - return; double t0 = PCU_Time(); EdgeReshaper es(a); ma::applyOperator(a,&es); @@ -837,17 +836,17 @@ void fixCrvElementShapes(Adapt* a) if ( ! count) break; prev_count = count; - //fixLargeAngles(a); // update this - //int numOpSuccess = fixLargeAngles(a); // update this */ - //PCU_Add_Ints(&numOpSuccess,1); - //if (PCU_Comm_Self() == 0) - // lion_oprint(1,"==> %d large angle fix operations succeeded.\n", numOpSuccess); - //markCrvBadQuality(a); + fixLargeAngles(a); // update this + /* int numOpSuccess = fixLargeAngles(a); // update this */ + /* PCU_Add_Ints(&numOpSuccess,1); */ + /* if (PCU_Comm_Self() == 0) */ + /* lion_oprint(1,"==> %d large angle fix operations succeeded.\n", numOpSuccess); */ + markCrvBadQuality(a); fixShortEdgeElements(a); // update this - //int numEdgeRemoved = fixShortEdgeElements(a); // update this */ - //PCU_Add_Ints(&numEdgeRemoved,1); - //if (PCU_Comm_Self() == 0) - // lion_oprint(1,"==> %d edges removal operations succeeded.\n", numEdgeRemoved); + /* int numEdgeRemoved = fixShortEdgeElements(a); // update this */ + /* PCU_Add_Ints(&numEdgeRemoved,1); */ + /* if (PCU_Comm_Self() == 0) */ + /* lion_oprint(1,"==> %d edges removal operations succeeded.\n", numEdgeRemoved); */ count = markCrvBadQuality(a); ++i; } while(count < prev_count && i < 6); // the second conditions is to make sure this does not take long diff --git a/ma/CMakeLists.txt b/ma/CMakeLists.txt index 6ef1bc036..68e72d2ea 100644 --- a/ma/CMakeLists.txt +++ b/ma/CMakeLists.txt @@ -56,8 +56,6 @@ set(HEADERS ma.h maInput.h maMesh.h - maMap.h - maAffine.h maSize.h maShape.h maTables.h diff --git a/ma/maShape.cc b/ma/maShape.cc index dbf42b3aa..458c9f4eb 100644 --- a/ma/maShape.cc +++ b/ma/maShape.cc @@ -814,7 +814,6 @@ void fixElementShapes(Adapt* a) return; double t0 = PCU_Time(); int count = markBadQuality(a); - printf(" number of bad shape inside fixElementShape %d \n", count); int originalCount = count; int prev_count; double time; diff --git a/ma/maSize.cc b/ma/maSize.cc index 64b64d199..c0fab50c7 100644 --- a/ma/maSize.cc +++ b/ma/maSize.cc @@ -429,11 +429,8 @@ struct LogAnisoSizeField : public MetricSizeField logMEval(f) { mesh = m; - // ADDED logMField = apf::createUserField(m, "ma_logM", apf::MATRIX, apf::getLagrange(1), &logMEval); - //logMField = apf::createUserField(m, "ma_logM", apf::MATRIX, - // apf::getLagrange(2), &logMEval); } ~LogAnisoSizeField() { @@ -442,9 +439,7 @@ struct LogAnisoSizeField : public MetricSizeField void init(Mesh* m, apf::Field* sizes, apf::Field* frames) { mesh = m; - //ADDED logMField = apf::createField(m, "ma_logM", apf::MATRIX, apf::getLagrange(1)); - //logMField = apf::createField(m, "ma_logM", apf::MATRIX, apf::getLagrange(2)); Entity* v; Iterator* it = m->begin(0); while ( (v = m->iterate(it)) ) { @@ -458,25 +453,6 @@ struct LogAnisoSizeField : public MetricSizeField 0 , 0 , s[2]); apf::setMatrix(logMField, v, 0, f * S * transpose(f)); } - m->end(it); - - //ADDED - /* - it = m->begin(1); - while ( (v = m->iterate(it)) ) { - Vector h; - Matrix f; - apf::getVector(sizes, v, 0, h); - apf::getMatrix(frames, v, 0, f); - Vector s(log(1/h[0]/h[0]), log(1/h[1]/h[1]), log(1/h[2]/h[2])); - Matrix S(s[0], 0 , 0, - 0 , s[1], 0, - 0 , 0 , s[2]); - apf::setMatrix(logMField, v, 0, f * S * transpose(f)); - } - m->end(it); - */ - } void getTransform( apf::MeshElement* me, From 3caaeb22bc652838ccb6fb29cc82b4a5dcdb3c95 Mon Sep 17 00:00:00 2001 From: Avinash Date: Wed, 2 Sep 2020 01:33:17 -0400 Subject: [PATCH 32/36] cleanup. --- crv/CMakeLists.txt | 1 - crv/crvAdapt.cc | 5 ++--- crv/crvShapeHandler.cc | 4 +--- ma/maRefine.cc | 2 +- 4 files changed, 4 insertions(+), 8 deletions(-) diff --git a/crv/CMakeLists.txt b/crv/CMakeLists.txt index 86f650ae0..302192c05 100644 --- a/crv/CMakeLists.txt +++ b/crv/CMakeLists.txt @@ -22,7 +22,6 @@ set(SOURCES crvTables.cc crvQuality.cc crvVtk.cc - crvBezierFields.cc ) # Package headers diff --git a/crv/crvAdapt.cc b/crv/crvAdapt.cc index 86214aabc..dffc2236b 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -17,8 +17,7 @@ #include #include #include -#include -#include + namespace crv { Adapt::Adapt(ma::Input* in) @@ -195,7 +194,7 @@ static void flagCleaner(crv::Adapt* a) } } -void getAllBezierFields(ma::Mesh* m, std::vector& fields) +static void getAllBezierFields(ma::Mesh* m, std::vector& fields) { for (int i = 0; i < m->countFields(); i++) { apf::FieldShape* fs = apf::getShape(m->getField(i)); diff --git a/crv/crvShapeHandler.cc b/crv/crvShapeHandler.cc index 7d4bf41ca..d8324bb74 100644 --- a/crv/crvShapeHandler.cc +++ b/crv/crvShapeHandler.cc @@ -19,7 +19,7 @@ #include #include #include -#include + namespace crv { static double measureLinearTriArea(ma::Mesh* m, ma::Entity* tri) @@ -207,8 +207,6 @@ class BezierTransfer : public ma::SolutionTransfer for (int k = 0; k < np; ++k) { point += nodes[k]*B(j+n-ni,k); } - //std::cout<<" shape ("<setPoint(newEntities[i],j,point); } } diff --git a/ma/maRefine.cc b/ma/maRefine.cc index 794cc5960..738fd5733 100644 --- a/ma/maRefine.cc +++ b/ma/maRefine.cc @@ -20,7 +20,7 @@ #include "maLayer.h" #include #include -#include + namespace ma { void addEdgePreAllocation(Refine* r, Entity* e, int counts[4]) From 03cd0a2d3142436bfac51e6f1a6641ecbd5cd3c7 Mon Sep 17 00:00:00 2001 From: Avinash Date: Tue, 8 Sep 2020 11:09:46 -0400 Subject: [PATCH 33/36] removed print blending order. --- crv/crvBezierFields.cc | 7 ------- 1 file changed, 7 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index c63a700b7..96bb546dc 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -86,13 +86,6 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) // revert back to at the end of this (so the rest of Bezier related // procedures work as intendedi) // - /* - int oldBlendings[apf::Mesh::TYPES]; - for (int i = 0; i < apf::Mesh::TYPES; i++) { - oldBlendings[i] = getBlendingOrder(i); - std::cout<<" blending order "<< oldBlendings[i] <<" for type "<< i< Date: Sun, 13 Sep 2020 22:21:35 -0400 Subject: [PATCH 34/36] cleanup. --- crv/crvBezierFields.cc | 8 -------- crv/crvBezierSolutionTransfer.cc | 9 --------- crv/crvShapeHandler.cc | 5 ++--- ma/maSolutionTransfer.cc | 16 ---------------- ma/maSolutionTransfer.h | 6 ++---- 5 files changed, 4 insertions(+), 40 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index 96bb546dc..bf50dc2d2 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -80,14 +80,6 @@ void convertInterpolationFieldPoints(apf::MeshEntity* e, void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) { - // TODO: to be completed - // Zero out the blending orders (since they are used internally - // in Bezier shapes (triangles and tets) but save the old ones to - // revert back to at the end of this (so the rest of Bezier related - // procedures work as intendedi) - // - //setBlendingOrder(apf::Mesh::TYPES, 0); - apf::FieldShape * fs = apf::getShape(f); int order = fs->getOrder(); diff --git a/crv/crvBezierSolutionTransfer.cc b/crv/crvBezierSolutionTransfer.cc index 110670887..053907468 100644 --- a/crv/crvBezierSolutionTransfer.cc +++ b/crv/crvBezierSolutionTransfer.cc @@ -184,10 +184,6 @@ class CrvBezierSolutionTransfer : public ma::SolutionTransfer return others.hasNodesOn(dimension); } - virtual const char* getTransferFieldName() { - return apf::getName(f); - } - virtual void onVertex( apf::MeshElement* parent, ma::Vector const& xi, @@ -414,11 +410,6 @@ ma::SolutionTransfer* setBezierSolutionTransfers( for (std::size_t i = 0; i < fields.size(); i++) { st->add(createBezierSolutionTransfer(fields[i], a)); } - std::vector trans = st->transfers; - for (std::size_t i = 0; i < trans.size(); i++) { - const char* name = trans[i]->getTransferFieldName(); - std::cout<<" field name added "<< name< A(n,np),B(n,np); getBezierTransformationMatrix(parentType,childType,P,A,vp); mth::multiply(Ai[apf::Mesh::typeDimension[childType]],A,B); - for (int j = 0; j < ni; ++j){ + for (int j = 0; j < ni; ++j) { apf::Vector3 point(0,0,0); - for (int k = 0; k < np; ++k) { + for (int k = 0; k < np; ++k) point += nodes[k]*B(j+n-ni,k); - } mesh->setPoint(newEntities[i],j,point); } } diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index 28a4aef98..dba3f1732 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -60,12 +60,6 @@ int SolutionTransfer::getTransferDimension() } return transferDimension; } - -const char* SolutionTransfer::getTransferFieldName() -{ - return 0; -} - FieldTransfer::FieldTransfer(apf::Field* f) { field = f; @@ -78,12 +72,6 @@ bool FieldTransfer::hasNodesOn(int dimension) { return shape->hasNodesIn(dimension); } - -const char* FieldTransfer::getTransferFieldName() -{ - return apf::getName(field); -} - void LinearTransfer::onVertex( apf::MeshElement* parent, Vector const& xi, @@ -288,10 +276,6 @@ class HighOrderTransfer : public SolutionTransfer { others.onCavity(oldElements,newEntities); } - const char* getTransferFieldName() - { - return verts.getTransferFieldName(); - } }; SolutionTransfer* createFieldTransfer(apf::Field* f) diff --git a/ma/maSolutionTransfer.h b/ma/maSolutionTransfer.h index 46df6fd84..15d183f59 100644 --- a/ma/maSolutionTransfer.h +++ b/ma/maSolutionTransfer.h @@ -77,7 +77,6 @@ class SolutionTransfer EntityArray& newEntities); /** \brief for internal MeshAdapt use */ int getTransferDimension(); - virtual const char* getTransferFieldName(); }; class FieldTransfer : public SolutionTransfer @@ -89,7 +88,6 @@ class FieldTransfer : public SolutionTransfer apf::Mesh* mesh; apf::FieldShape* shape; apf::NewArray value; - virtual const char* getTransferFieldName(); }; class LinearTransfer : public FieldTransfer @@ -156,7 +154,7 @@ class SolutionTransfers : public SolutionTransfer virtual bool hasNodesOn(int dimension); virtual void onVertex( apf::MeshElement* parent, - Vector const& xi, + Vector const& xi, Entity* vert); virtual void onRefine( Entity* parent, @@ -164,7 +162,7 @@ class SolutionTransfers : public SolutionTransfer virtual void onCavity( EntityArray& oldElements, EntityArray& newEntities); - //private + private: typedef std::vector Transfers; Transfers transfers; }; From cfe92b423180fba64ab11a3367707bf2b0d0f973 Mon Sep 17 00:00:00 2001 From: Avinash Date: Tue, 15 Sep 2020 04:26:30 -0400 Subject: [PATCH 35/36] regression test for bezier field. --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index af2bc7936..e84a30c06 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -184,6 +184,7 @@ test_exe_func(create_mis create_mis.cc) test_exe_func(fieldReduce fieldReduce.cc) test_exe_func(test_integrator test_integrator.cc) test_exe_func(test_matrix_gradient test_matrix_grad.cc) +test_exe_func(bezierField bezierField.cc) if(ENABLE_DSP) test_exe_func(graphdist graphdist.cc) test_exe_func(moving moving.cc) From e73749b0f60f4d813ad48f63348f0adcbe487a61 Mon Sep 17 00:00:00 2001 From: Avinash Date: Thu, 24 Sep 2020 14:08:29 -0400 Subject: [PATCH 36/36] convert interpolating to bezier fix for blended entities. --- crv/crvBezierFields.cc | 5 +++-- crv/crvCurveMesh.cc | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc index bf50dc2d2..cbf488a96 100644 --- a/crv/crvBezierFields.cc +++ b/crv/crvBezierFields.cc @@ -95,7 +95,6 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) int startDim = md - (blendingOrder > 0); for(int d = startDim; d >= 1; --d){ - //for(int d = md; d >= 1; d--){ if(!fs->hasNodesIn(d)) continue; int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); @@ -112,8 +111,10 @@ void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f) } for( int d = 2; d <= md; ++d){ + std::cout<<" blendingorder of dimension "<hasNodesIn(d) || - getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; + !getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); apf::NewArray c; diff --git a/crv/crvCurveMesh.cc b/crv/crvCurveMesh.cc index 5cc39bbfb..9c31ef418 100644 --- a/crv/crvCurveMesh.cc +++ b/crv/crvCurveMesh.cc @@ -152,7 +152,7 @@ void BezierCurver::convertInterpolatingToBezier() // triangles and tetrahedra for(int d = 2; d <= md; ++d){ if(!fs->hasNodesIn(d) || - getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; + !getBlendingOrder(apf::Mesh::simplexTypes[d])) continue; int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes(); int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]); apf::NewArray c;