diff --git a/crv/CMakeLists.txt b/crv/CMakeLists.txt index 992589d23..302192c05 100644 --- a/crv/CMakeLists.txt +++ b/crv/CMakeLists.txt @@ -6,6 +6,8 @@ set(SOURCES crvBezier.cc crvBezierPoints.cc crvBezierShapes.cc + crvBezierFields.cc + crvBezierSolutionTransfer.cc crvBlended.cc crvCurveMesh.cc crvElevation.cc diff --git a/crv/crv.h b/crv/crv.h index 60e8c2c1a..12a8d4850 100644 --- a/crv/crv.h +++ b/crv/crv.h @@ -11,6 +11,7 @@ #include "apfMesh2.h" #include "apfShape.h" #include +#include #include #include #include @@ -22,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; @@ -37,6 +41,13 @@ int getBlendingOrder(const int type); /** \brief count invalid elements of the mesh */ int countNumberInvalidElements(apf::Mesh2* m); +/** \ 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 \details P is the order, S is the space dimension, different from the mesh dimension, used to distinguish between planar 2D @@ -203,6 +214,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, 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 3b6baec15..dffc2236b 100644 --- a/crv/crvAdapt.cc +++ b/crv/crvAdapt.cc @@ -194,6 +194,16 @@ static void flagCleaner(crv::Adapt* a) } } +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)); + std::string name = fs->getName(); + if (name == std::string("Bezier")) + fields.push_back(m->getField(i)); + } +} + void adapt(ma::Input* in) { std::string name = in->mesh->getShape()->getName(); @@ -205,6 +215,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::setBezierSolutionTransfers(allFields, a); + a->solutionTransfer = in->solutionTransfer; + ma::preBalance(a); fixInvalidElements(a); diff --git a/crv/crvBezierFields.cc b/crv/crvBezierFields.cc new file mode 100644 index 000000000..cbf488a96 --- /dev/null +++ b/crv/crvBezierFields.cc @@ -0,0 +1,136 @@ +/* + * 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 "crvBezierShapes.h" +#include "crvMath.h" +#include "crvShape.h" +#include "crvTables.h" +#include "crvQuality.h" +#include +#include +#include +#include +#include +#include + +namespace crv { + +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::MeshElement* me = + apf::createMeshElement(f->getMesh(),e); + apf::Element* elem = + apf::createElement(f,me); + + 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); + apf::destroyMeshElement(me); +} + +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]); + //blendingOrder = 0; + // 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); + } + + for( int d = 2; d <= md; ++d){ + std::cout<<" blendingorder of dimension "<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 new file mode 100644 index 000000000..053907468 --- /dev/null +++ b/crv/crvBezierSolutionTransfer.cc @@ -0,0 +1,416 @@ +/* + 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 "crv.h" +#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 +#include +#include + +#include +#include +namespace crv { + +static void setLinearEdgeFieldPoints(ma::Mesh* m, apf::Field* f, + ma::Entity* edge) +{ + apf::Vector3 xi,points[2]; + apf::MeshEntity* verts[2]; + apf::FieldShape* fs = apf::getShape(f); + + int ni = fs->countNodesOn(apf::Mesh::EDGE); + m->getDownward(edge,0,verts); + 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])); + } +} + +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; + } + } + } + } +} + +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; + + 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); + apf::getVector(e, xin, xyz); + fval = (xyz-point); + + xinew = xin - transpose(Jinv)*fval; + + 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 xin = xinew; + } + else { + xin = xinew; + } + } + apf::destroyElement(e); + } + + return (elemNum); +} + +class CrvBezierSolutionTransfer : public ma::SolutionTransfer +{ + public: + Adapt* adapt; + ma::CavityTransfer others; + + CrvBezierSolutionTransfer(apf::Field* field, Adapt* a): + adapt(a),others(field) + { + 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]); + } + } + virtual bool hasNodesOn(int dimension) + { + return others.hasNodesOn(dimension); + } + + 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; + + mesh->getDownward(parent, 0, parentVerts); + mesh->getDownward(parent, 1, parentEdges); + + int ne = apf::Mesh::adjacentCount[parentType][1]; + + 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; + if ( ma::getFlag(adapt,parentEdges[i], ma::BAD_QUALITY) ) { + useLinear = true; + } + } + + 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 = 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 { + + 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],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); + } + } + } + } + } + } + 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::NewArray elems(n); + apf::NewArray elemsF(n); + for (int i = 0; i < n; i++) { + elems[i] = apf::createMeshElement(mesh, oldElements[i]); + } + + entNum = getBestElement(n, mesh, + &elems[0], xyz, xi); + + PCU_ALWAYS_ASSERT(entNum >= 0); + + apf::MeshElement* meshElemP = apf::createMeshElement( + mesh, oldElements[entNum]); + apf::Vector3 xxyy, xxzz, val2; + apf::mapLocalToGlobal(meshElemP, xi, xxyy); + + //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); + } + + 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]); + + if (!ne) continue; + + apf::MeshElement* me = apf::createMeshElement(mesh, newEntities[i]); + apf::Element* e = apf::createElement(mesh->getCoordinateField(), me); + apf::Vector3 point, xinew; + + 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; + + 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); + } + } + } + + 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)); + } + return st; +} + +} 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; diff --git a/crv/crvShapeHandler.cc b/crv/crvShapeHandler.cc index 9ec487ec0..e63cbcad1 100644 --- a/crv/crvShapeHandler.cc +++ b/crv/crvShapeHandler.cc @@ -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); @@ -198,10 +199,10 @@ class BezierTransfer : public ma::SolutionTransfer mesh->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){ + for (int j = 0; j < ni; ++j) { apf::Vector3 point(0,0,0); for (int k = 0; k < np; ++k) point += nodes[k]*B(j+n-ni,k); diff --git a/ma/maRefine.cc b/ma/maRefine.cc index 4325186fb..738fd5733 100644 --- a/ma/maRefine.cc +++ b/ma/maRefine.cc @@ -345,14 +345,16 @@ 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) - 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) + 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]); + } } void forgetNewEntities(Refine* r) diff --git a/ma/maSolutionTransfer.cc b/ma/maSolutionTransfer.cc index 89b564e82..dba3f1732 100644 --- a/ma/maSolutionTransfer.cc +++ b/ma/maSolutionTransfer.cc @@ -10,10 +10,8 @@ #include "maSolutionTransfer.h" #include "maAffine.h" #include "maMap.h" -#include -#include -#include - +#include "apfField.h" +#include namespace ma { SolutionTransfer::~SolutionTransfer() @@ -62,143 +60,177 @@ 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); +} +void LinearTransfer::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); +} -class CavityTransfer : public FieldTransfer +CavityTransfer::CavityTransfer(apf::Field* f): + FieldTransfer(f) { - public: - 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) { - minDim = getMinimumDimension(getShape(f)); + bestValue = value; + bestI = i; + bestXi = xi; } - void transferToNodeIn( - apf::Element* elem, - apf::Node const& node, - Vector const& elemXi) + } + 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) { - apf::getComponents(elem,elemXi,&(value[0])); - apf::setComponents(field,node.entity,node.node,&(value[0])); + apf::Node node(newEntities[i],j); + transferToNode(n,&(elems[0]),&(elemInvMaps[0]),node); } - 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; - } + } + for (int i = 0; i < n; ++i) + apf::destroyElement(elems[i]); +} +void CavityTransfer::onRefine( + Entity* parent, + EntityArray& 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); } - 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 CavityTransfer::onCavity( + EntityArray& oldElements, + EntityArray& newEntities) +{ + apf::FieldShape *fs = apf::getShape(field); + std::string name = fs->getName(); + bool fConstant = false; + if (name == std::string("Constant_3")) + fConstant = true; + + double sumOld = 0.; + if (fConstant) { + 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; } - 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); - } + } + transfer(oldElements.getSize(),&(oldElements[0]),newEntities); + + double sumNew = 0.; + double sumNewVal = 0.; + if (fConstant) { + 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]); + sumNew += val*apf::measure(mesh, newEntities[i]); + //sumNew += val; } - for (int i = 0; i < n; ++i) - apf::destroyElement(elems[i]); } - virtual void onRefine( - Entity* parent, - EntityArray& newEntities) - { - transfer(1,&parent,newEntities); + 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); + if (sumNew == 0) continue; + apf::setScalar(field, newEntities[i], 0, val*sumOld/sumNew); + } } - virtual void onCavity( - EntityArray& oldElements, - EntityArray& newEntities) - { - transfer(oldElements.getSize(),&(oldElements[0]),newEntities); + 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]); + } } - private: - int minDim; -}; + } +} /* hmm... could use multiple inheritance here, but that creates a "diamond problem": @@ -313,7 +345,9 @@ AutoSolutionTransfer::AutoSolutionTransfer(Mesh* m) for (int i = 0; i < m->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 a80fe90e5..15d183f59 100644 --- a/ma/maSolutionTransfer.h +++ b/ma/maSolutionTransfer.h @@ -14,10 +14,16 @@ \brief MeshAdapt solution transfer interface */ #include +#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 +79,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 @@ -92,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, 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)