From a1f86f7b3bde90d7ffae5bc1521cb75eaeb2c15a Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 30 Jan 2019 21:07:22 +0100 Subject: [PATCH 01/17] [temp] integration 1d3d --- .../embedded/couplingmanager1d3d.hh | 313 +++++++++++------- 1 file changed, 187 insertions(+), 126 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 53e647e594..5163cbec3a 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -26,6 +26,7 @@ #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_HH #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_HH +#include #include #include @@ -39,6 +40,7 @@ #include #include #include +#include namespace Dumux { @@ -972,89 +974,101 @@ public: // intersect the bounding box trees this->glueGrids(); + // reserve memory for data this->pointSourceData().reserve(this->glue().size()); this->averageDistanceToBulkCell().reserve(this->glue().size()); - const Scalar kernelWidth = getParam("MixedDimension.KernelWidth"); + + // generate a bunch of random vectors and values for + // Monte-carlo integration on the cylinder defined by line and radius + const int samples = getParam("MixedDimension.NumMCSamples"); + CylinderIntegration cylIntegration(samples); + static const auto kernelWidthFactor = getParam("MixedDimension.KernelWidthFactor"); + for (const auto& is : intersections(this->glue())) { // all inside elements are identical... const auto& inside = is.targetEntity(0); // get the intersection geometry for integrating over it const auto intersectionGeometry = is.geometry(); - // get the Gaussian quadrature rule for the local intersection - const auto& quad = Dune::QuadratureRules::rule(intersectionGeometry.type(), order); - const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside); + const auto lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(inside); - // iterate over all quadrature points and place a source - // for 1d: make a new point source - // for 3d: make a new kernel volume source - for (auto&& qp : quad) + // for each intersection integrate kernel and add: + // * 1d: a new point source + // * 3d: a new kernel volume source + const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); + const auto kernelWidth = kernelWidthFactor*radius; + + + // we can have multiple 3d elements per 1d intersection, for evaluation taking one of them should be fine + for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) { - // compute the coupling stencils - for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) - { - const auto& outside = is.domainEntity(outsideIdx); - const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside); + // compute source id + const auto id = this->idCounter_++; - // each quadrature point will be a point source for the sub problem - const auto globalPos = intersectionGeometry.global(qp.position()); - const auto id = this->idCounter_++; - const auto qpweight = qp.weight(); - const auto ie = intersectionGeometry.integrationElement(qp.position()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, std::vector({lowDimElementIdx})); - this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); - computeBulkSource(globalPos, kernelWidth, id, lowDimElementIdx, bulkElementIdx, qpweight*ie/is.numDomainNeighbors()); + // place a point source at the intersection center + const auto center = intersectionGeometry.center(); + this->pointSources(lowDimIdx).emplace_back(center, id, /*weight=*/1.0, intersectionGeometry.volume(), std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); - // pre compute additional data used for the evaluation of - // the actual solution dependent source term - PointSourceData psData; + const auto& outside = is.domainEntity(outsideIdx); + const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside); - if constexpr (isBox()) - { - using ShapeValues = std::vector >; - const auto lowDimGeometry = lowDimGridGeometry.element(lowDimElementIdx).geometry(); - ShapeValues shapeValues; - this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues); - psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); - } - else - { - psData.addLowDimInterpolation(lowDimElementIdx); - } + // compute the weights for the bulk volume sources + static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); + if (useCylinderIntegration) + cylIntegration.setGeometry(intersectionGeometry.corner(0), intersectionGeometry.corner(1), kernelWidth); + computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); - // add data needed to compute integral over the circle - if constexpr (isBox()) - { - using ShapeValues = std::vector >; - const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); - ShapeValues shapeValues; - this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, globalPos, shapeValues); - psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); - } - else - { - psData.addBulkInterpolation(bulkElementIdx); - } + // pre compute additional data used for the evaluation of + // the actual solution dependent source term + PointSourceData psData; - // publish point source data in the global vector - this->pointSourceData().emplace_back(std::move(psData)); + // lowdim interpolation (evaluate at center) + if (isBox()) + { + using ShapeValues = std::vector >; + const auto lowDimGeometry = this->problem(lowDimIdx).fvGridGeometry().element(lowDimElementIdx).geometry(); + ShapeValues shapeValues; + this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).fvGridGeometry(), lowDimGeometry, center, shapeValues); + psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); + } + else + { + psData.addLowDimInterpolation(lowDimElementIdx); + } - // compute average distance to bulk cell - this->averageDistanceToBulkCell().push_back(averageDistancePointGeometry(globalPos, outside.geometry())); + // bulk interpolation (evaluate at center) + if (isBox()) + { + using ShapeValues = std::vector >; + const auto bulkGeometry = this->problem(bulkIdx).fvGridGeometry().element(bulkElementIdx).geometry(); + ShapeValues shapeValues; + this->getShapeValues(bulkIdx, this->problem(bulkIdx).fvGridGeometry(), bulkGeometry, center, shapeValues); + psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); + } + else + { + psData.addBulkInterpolation(bulkElementIdx); + } - // export the lowdim coupling stencil - // we insert all vertices / elements and make it unique later - if (isBox()) - { - const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); - this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), - vertices.begin(), vertices.end()); - } - else - { - this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); - } + // publish point source data in the global vector + this->pointSourceData().emplace_back(std::move(psData)); + + // compute average distance to bulk cell + this->averageDistanceToBulkCell().push_back(this->computeDistance(outside.geometry(), center)); + + // export the lowdim coupling stencil + // we insert all vertices / elements and make it unique later + if (isBox()) + { + const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); + this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), + vertices.begin(), vertices.end()); + } + else + { + this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); } } } @@ -1166,85 +1180,132 @@ public: // \} private: - //! TODO: How to optimize this? - void computeBulkSource(const GlobalPosition& globalPos, const Scalar kernelWidth, - std::size_t id, GridIndex lowDimElementIdx, GridIndex coupledBulkElementIdx, - Scalar pointSourceWeight) + + template + void computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, + std::size_t id, std::size_t lowDimElementIdx, std::size_t coupledBulkElementIdx, + const CylIntegration& cylIntegration, std::size_t embeddings) { - // make sure it is mass conservative - // i.e. the point source in the 1d domain needs to have the exact same integral as the distributed - // kernel source integral in the 3d domain. Correct the integration formula by balancing the error - // by scaling the kernel with the checksum inverse - Scalar checkSum = 0.0; - std::vector mask(this->gridView(bulkIdx).size(0), false); - for (const auto& element : elements(this->gridView(bulkIdx))) - { - Scalar weight = 0.0; - const auto geometry = element.geometry(); - const auto& quad = Dune::QuadratureRules::rule(geometry.type(), 3); - for (auto&& qp : quad) - { - const auto qpweight = qp.weight(); - const auto ie = geometry.integrationElement(qp.position()); - weight += evalKernel(globalPos, geometry.global(qp.position()), kernelWidth)*qpweight*ie; - } + // Monte-carlo integration on the cylinder defined by line and radius + const auto numBulkElements = this->gridView(bulkIdx).size(0); + std::vector weights(numBulkElements, 0.0); + std::vector visited(numBulkElements, false); - if (weight > 1e-13) - { - const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); - const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(element); - bulkSourceIds_[bulkElementIdx].push_back(id); - bulkSourceWeights_[bulkElementIdx].push_back(weight*pointSourceWeight); - mask[bulkElementIdx] = true; + Scalar integral = 0.0; - // add lowDim dofs that the source is related to to the bulk stencil - if (isBox()) + static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); + if (useCylinderIntegration) + { + for (int i = 0; i < cylIntegration.size(); ++i) + { + const auto point = cylIntegration.getIntegrationPoint(i); + const auto weight = evalKernel(line.corner(0), line.corner(1), point, radius, kernelWidth)*cylIntegration.integrationElement(i); + integral += weight; + + static const auto min = getParam("Tissue.Grid.LowerLeft"); + static const auto max = getParam("Tissue.Grid.UpperRight"); + static const auto cells = getParam("Tissue.Grid.Cells"); + const auto is = intersectingEntityCartesian(point, min, max, cells); + // const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); + if (is.first) { - const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx); - this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), - vertices.begin(), vertices.end()); - + const auto bulkIdx = is.second; + visited[bulkIdx] = true; + weights[bulkIdx] += weight; } - else + } + } + // integrate over 3d domain using standard quadrature rules + else + { + static const int maxLevel = getParam("MixedDimension.MaxOctreeLevel"); + CylinderDomain cylinder(line.corner(0), line.corner(1), kernelWidth); + const auto kernelFunc = [&](const auto& p){ return evalKernel(line.corner(0), line.corner(1), p, radius, kernelWidth); }; + for (const auto& element : elements(this->gridView(bulkIdx))) + { + const auto bulkElementIdx = this->problem(bulkIdx).fvGridGeometry().elementMapper().index(element); + const auto geometry = element.geometry(); + OctreeIntegrator integrator(geometry.corner(0), geometry.corner(7), maxLevel); + const auto weight = integrator.integrate(kernelFunc, cylinder); + if (weight > 1e-35) { - this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); + visited[bulkElementIdx] = true; + weights[bulkElementIdx] += weight; + integral += weight; } - - // tpfa - extendedSourceStencil_.stencil()[bulkElementIdx].push_back(coupledBulkElementIdx); - - // compute check sum -> should sum up to 1.0 to be mass conservative - checkSum += weight; } } - for (GridIndex eIdx = 0; eIdx < bulkSourceWeights_.size(); ++eIdx) - if (mask[eIdx]) bulkSourceWeights_[eIdx].back() /= checkSum; + const auto length = line.volume(); + std::cout << "integration error in %: " << std::abs(length-integral)/length*100 << std::endl; + // const auto correctionFactor = length/integral; + for (auto bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) + { + if (!visited[bulkElementIdx]) + continue; + + bulkSourceIds_[bulkElementIdx].push_back(id); + bulkSourceWeights_[bulkElementIdx].push_back(weights[bulkElementIdx]/embeddings);//*correctionFactor); + + if (isBox()) + { + const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx); + this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), + vertices.begin(), vertices.end()); + + } + else + { + this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); + } - // balance error of the quadrature rule -> TODO: what to do at boundaries - // const auto diff = 1.0 - checkSum; - // std::cout << "Integrated kernel with integration error of " << diff << std::endl; + // the extended source stencil, every 3d element with a source is coupled to the element where the 3d quantities are measured + extendedSourceStencil_.stencil()[bulkElementIdx].push_back(coupledBulkElementIdx); + } } - //! an isotropic cubic kernel with derivatives 0 at r=origin and r=width and domain integral 1 - inline Scalar evalKernel(const GlobalPosition& origin, - const GlobalPosition& pos, - const Scalar width) const noexcept + //! a cylindrical kernel around the segment a->b + inline Scalar evalKernel(const GlobalPosition& a, + const GlobalPosition& b, + const GlobalPosition& point, + const Scalar R, + const Scalar rho) const noexcept { - const auto r = (pos-origin).two_norm(); - const auto r2 = r*r; - const auto r3 = r2*r; + // projection of point onto line a + t*(b-a) + const auto ab = b - a; + const auto t = (point - a)*ab/ab.two_norm2(); - if (r > width) + // return 0 if we are outside cylinder + if (t < 0.0 || t > 1.0) return 0.0; - const Scalar w2 = width*width; - const Scalar w3 = w2*width; - const Scalar k = 15.0/(4*M_PI*w3); - const Scalar a = 2.0/w3; - const Scalar b = 3.0/w2; + // compute distance + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - point).two_norm(); + + if (r > rho) + return 0.0; - return k*(a*r3 - b*r2 + 1.0); + // custom kernel (constant) + return 1.0/(M_PI*rho*rho); + + // custom kernel (linear) + // const auto lnR = std::log(R); + // const auto lnRho = std::log(rho); + // return (18.0*r*lnR - 18.0*r*lnRho + 9*r - 12.0*rho*lnR + 12.0*rho*lnRho - 4*rho)/(2*M_PI*rho*rho*rho); + + // custom kernel (quadratic) + // const auto lnR = std::log(R); + // const auto lnRho = std::log(rho); + // return (8.0*r*r*lnR - 8.0*r*r*lnRho + 4.0*r*r - 4.0*rho*rho*lnR + 4.0*rho*rho*lnRho - rho*rho)/(M_PI*rho*rho*rho*rho); + + // cubic kernel + // const auto r2 = r*r; + // const auto r3 = r2*r; + // const Scalar w2 = rho*rho; + // const Scalar w3 = w2*rho; + // const Scalar k = 10.0/3.0/M_PI/w2; + // return k*(2.0*r3/w3 - 3.0*r2/w2 + 1.0); } //! the extended source stencil object for kernel coupling -- GitLab From 448baa79c78e25248624dc225aaac01f32c3c8c6 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 12 Feb 2019 16:27:06 +0100 Subject: [PATCH 02/17] [temp] Couplingmanager 1d3d --- .../embedded/couplingmanager1d3d.hh | 192 +++++++++++++----- 1 file changed, 137 insertions(+), 55 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 5163cbec3a..69cdb13a3c 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -30,6 +30,7 @@ #include #include +#include #include #include @@ -977,11 +978,13 @@ public: // reserve memory for data this->pointSourceData().reserve(this->glue().size()); this->averageDistanceToBulkCell().reserve(this->glue().size()); + fluxScalingFactor_.reserve(this->glue().size()); // generate a bunch of random vectors and values for // Monte-carlo integration on the cylinder defined by line and radius - const int samples = getParam("MixedDimension.NumMCSamples"); - CylinderIntegration cylIntegration(samples); + static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); + const int samples = useCylinderIntegration ? getParam("MixedDimension.NumMCSamples") : 1; + CylinderIntegration cylIntegration(samples); static const auto kernelWidthFactor = getParam("MixedDimension.KernelWidthFactor"); for (const auto& is : intersections(this->glue())) @@ -999,11 +1002,18 @@ public: const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); const auto kernelWidth = kernelWidthFactor*radius; + const auto a = intersectionGeometry.corner(0); + const auto b = intersectionGeometry.corner(1); // we can have multiple 3d elements per 1d intersection, for evaluation taking one of them should be fine for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) { // compute source id + // for each id we have + // (1) a point source for the 1d domain + // (2) a bulk source weight for the each element in the 3d domain (the fraction of the total source/sink for the 1d element with the point source) + // TODO: i.e. one integration of the kernel should be enough (for each of the weights it's weight/embeddings) + // (3) the flux scaling factor for each outside element, i.e. each id const auto id = this->idCounter_++; // place a point source at the intersection center @@ -1015,9 +1025,8 @@ public: const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside); // compute the weights for the bulk volume sources - static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); if (useCylinderIntegration) - cylIntegration.setGeometry(intersectionGeometry.corner(0), intersectionGeometry.corner(1), kernelWidth); + cylIntegration.setGeometry(a, b, kernelWidth); computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); // pre compute additional data used for the evaluation of @@ -1055,8 +1064,22 @@ public: // publish point source data in the global vector this->pointSourceData().emplace_back(std::move(psData)); - // compute average distance to bulk cell - this->averageDistanceToBulkCell().push_back(this->computeDistance(outside.geometry(), center)); + // compute the average flux scaling factor + const auto outsideGeometry = outside.geometry(); + static const int fluxScalingFactorIntOrder = getParam("MixedDimension.FluxScalingIntegrationOrder", 5); + const auto& quadBulk = Dune::QuadratureRules::rule(outsideGeometry.type(), fluxScalingFactorIntOrder); + // iterate over all quadrature points + Scalar avgMinDist = 0.0; + for (auto&& qpBulk : quadBulk) + { + const auto globalPosBulk = outsideGeometry.global(qpBulk.position()); + const auto d = computeMinDistance_(a, b, globalPosBulk); + avgMinDist += d*qpBulk.weight(); + } + + // add it to the internal data vector + this->averageDistanceToBulkCell().push_back(avgMinDist); + fluxScalingFactor_.push_back(this->problem(bulkIdx).fluxScalingFactor(avgMinDist, radius)); // export the lowdim coupling stencil // we insert all vertices / elements and make it unique later @@ -1170,17 +1193,21 @@ public: } //! return all source ids for a bulk elements - const std::vector bulkSourceIds(GridIndex eIdx) const - { return bulkSourceIds_[eIdx]; } + const std::vector& bulkSourceIds(GridIndex eIdx, std::size_t scvIdx = 0) const + { return bulkSourceIds_[eIdx][scvIdx]; } //! return all source ids for a bulk elements - const std::vector bulkSourceWeights(GridIndex eIdx) const - { return bulkSourceWeights_[eIdx]; } + const std::vector& bulkSourceWeights(GridIndex, std::size_t scvIdx = 0) const + { return bulkSourceWeights_[eIdx][scvIdx]; } + + //! The flux scaling factor for a source with id + Scalar fluxScalingFactor(std::size_t id) const + { return fluxScalingFactor_[id]; } // \} private: - + //! compute the kernel distributed sources and add stencils template void computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, std::size_t id, std::size_t lowDimElementIdx, std::size_t coupledBulkElementIdx, @@ -1188,11 +1215,11 @@ private: { // Monte-carlo integration on the cylinder defined by line and radius const auto numBulkElements = this->gridView(bulkIdx).size(0); - std::vector weights(numBulkElements, 0.0); - std::vector visited(numBulkElements, false); + constexpr std::size_t numScv = isBox() ? 1<> weights(numBulkElements, std::array{}); + std::vector> visited(numBulkElements, std::bitset{}); Scalar integral = 0.0; - static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); if (useCylinderIntegration) { @@ -1209,9 +1236,24 @@ private: // const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); if (is.first) { - const auto bulkIdx = is.second; - visited[bulkIdx] = true; - weights[bulkIdx] += weight; + const auto bulkElementIdx = is.second; + if (isBox()) + { + const auto dist = point-min; + const auto diag = max-min; + std::bitset<3> corner; + for (int i = 0; i < 3; ++i) + corner[i] = bool(int(std::round(dist[i]/diag[i]))); + + std::size_t scvIdx = corner.to_ulong(); + visited[bulkElementIdx][scvIdx] = true; + weights[bulkElementIdx][scvIdx] += weight; + } + else + { + visited[bulkElementIdx][0] = true; + weights[bulkElementIdx][0] += weight; + } } } } @@ -1229,37 +1271,74 @@ private: const auto weight = integrator.integrate(kernelFunc, cylinder); if (weight > 1e-35) { - visited[bulkElementIdx] = true; - weights[bulkElementIdx] += weight; + visited[bulkElementIdx][0] = true; + weights[bulkElementIdx][0] += weight; integral += weight; } } } - const auto length = line.volume(); - std::cout << "integration error in %: " << std::abs(length-integral)/length*100 << std::endl; + // const auto length = line.volume(); + // std::cout << "integration error in %: " << std::abs(length-integral)/length*100 << std::endl; // const auto correctionFactor = length/integral; - for (auto bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) + for (const auto& element : elements(this->problem(bulkIdx).fvGridGeometry().gridView())) { - if (!visited[bulkElementIdx]) + const auto bulkElementIdx = this->problem(bulkIdx).fvGridGeometry().elementMapper().index(element); + if (!isBox() && !visited[bulkElementIdx][0]) continue; - bulkSourceIds_[bulkElementIdx].push_back(id); - bulkSourceWeights_[bulkElementIdx].push_back(weights[bulkElementIdx]/embeddings);//*correctionFactor); - - if (isBox()) + if (isBox()) { - const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx); - this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), - vertices.begin(), vertices.end()); + auto fvGeometry = localView(this->problem(bulkIdx).fvGridGeometry()); + fvGeometry.bindElement(element); + for (const auto& scv : scvs(fvGeometry)) + { + const auto scvIdx = scv.indexInElement(); + if (!visited[bulkElementIdx][scvIdx]) + continue; + + bulkSourceIds_[bulkElementIdx][scvIdx].push_back(id); + bulkSourceWeights_[bulkElementIdx][scvIdx].push_back(weights[bulkElementIdx][scvIdx]/embeddings);//*correctionFactor); + addBulkSourceStencils_(bulkElementIdx, lowDimElementIdx, coupledBulkElementIdx); + } } else { - this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); + bulkSourceIds_[bulkElementIdx][0].push_back(id); + bulkSourceWeights_[bulkElementIdx][0].push_back(weights[bulkElementIdx][0]/embeddings);//*correctionFactor); + addBulkSourceStencils_(bulkElementIdx, lowDimElementIdx, coupledBulkElementIdx); } + } + + } + + //! add additional stencil entries for the bulk element + void addBulkSourceStencils_(std::size_t bulkElementIdx, std::size_t coupledLowDimElementIdx, std::size_t coupledBulkElementIdx) + { + // add the lowdim element to the coupling stencil of this bulk element + if (isBox()) + { + const auto& vertices = this->vertexIndices(lowDimIdx, coupledLowDimElementIdx); + this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), + vertices.begin(), vertices.end()); - // the extended source stencil, every 3d element with a source is coupled to the element where the 3d quantities are measured + } + else + { + this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(coupledLowDimElementIdx); + } + + // the extended source stencil, every 3d element with a source is coupled to + // the element/dofs where the 3d quantities are measured + if (isBox()) + { + const auto& vertices = this->vertexIndices(bulkIdx, coupledBulkElementIdx); + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + vertices.begin(), vertices.end()); + } + else + { extendedSourceStencil_.stencil()[bulkElementIdx].push_back(coupledBulkElementIdx); } } @@ -1269,7 +1348,7 @@ private: const GlobalPosition& b, const GlobalPosition& point, const Scalar R, - const Scalar rho) const noexcept + const Scalar rho) const { // projection of point onto line a + t*(b-a) const auto ab = b - a; @@ -1286,26 +1365,27 @@ private: if (r > rho) return 0.0; - // custom kernel (constant) - return 1.0/(M_PI*rho*rho); - - // custom kernel (linear) - // const auto lnR = std::log(R); - // const auto lnRho = std::log(rho); - // return (18.0*r*lnR - 18.0*r*lnRho + 9*r - 12.0*rho*lnR + 12.0*rho*lnRho - 4*rho)/(2*M_PI*rho*rho*rho); - - // custom kernel (quadratic) - // const auto lnR = std::log(R); - // const auto lnRho = std::log(rho); - // return (8.0*r*r*lnR - 8.0*r*r*lnRho + 4.0*r*r - 4.0*rho*rho*lnR + 4.0*rho*rho*lnRho - rho*rho)/(M_PI*rho*rho*rho*rho); - - // cubic kernel - // const auto r2 = r*r; - // const auto r3 = r2*r; - // const Scalar w2 = rho*rho; - // const Scalar w3 = w2*rho; - // const Scalar k = 10.0/3.0/M_PI/w2; - // return k*(2.0*r3/w3 - 3.0*r2/w2 + 1.0); + static const std::string type = getParam("MixedDimension.KernelType"); + if (type == "Const") + return 1.0/(M_PI*rho*rho); + else if (type == "Cubic") + { + const auto r2 = r*r; const auto r3 = r*r*r; + const auto rho2 = rho*rho; const auto rho3 = rho2*rho; + return 10.0*(2.0*r3/rho3 - 3.0*r2/rho2 + 1.0)/(3.0*M_PI*rho2); + } + + DUNE_THROW(Dune::NotImplemented, "Kernel type: " << type); + } + + inline Scalar computeMinDistance_(const GlobalPosition& a, const GlobalPosition& b, const GlobalPosition& x3d) const + { + // r is the minimum distance to the line through a and b + const auto ab = b - a; + const auto t = (x3d - a)*ab/ab.two_norm2(); + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - x3d).two_norm(); + return r; } //! the extended source stencil object for kernel coupling @@ -1313,9 +1393,11 @@ private: //! vector for the volume fraction of the lowdim domain in the bulk domain cells std::vector lowDimVolumeInBulkElement_; //! kernel sources to integrate for each bulk element - std::vector> bulkSourceIds_; + std::vector, isBox() ? 1<> bulkSourceIds_; //! the integral of the kernel for each point source / integration point, i.e. weight for the source - std::vector> bulkSourceWeights_; + std::vector, isBox() ? 1<> bulkSourceWeights_; + //! the flux scaling factor for the respective source id + std::vector fluxScalingFactor_; }; } // end namespace Dumux -- GitLab From aff5af959725122da42c589c5345144b070ac901 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 12 Feb 2019 18:10:03 +0100 Subject: [PATCH 03/17] [embedded] Improve average robustness by considering facet/edge cuts on the circle average --- .../embedded/couplingmanager1d3d.hh | 31 +++++++++---------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 69cdb13a3c..fc4581ad85 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -336,11 +336,12 @@ public: const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx); const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0); const auto circleAvgWeight = 2*M_PI*radius/numIp; + const auto integrationElement = lowDimGeometry.integrationElement(qp.position()); + const auto qpweight = qp.weight(); const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); std::vector> circleStencil; circleStencil.reserve(circlePoints.size()); - // for box std::vector>*> circleCornerIndices; using ShapeValues = std::vector >; @@ -354,24 +355,24 @@ public: if (circleBulkElementIndices.empty()) continue; - ++insideCirclePoints; const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices.size(); for (const auto bulkElementIdx : circleBulkElementIndices) { circleStencil.push_back(bulkElementIdx); circleIpWeight.push_back(localCircleAvgWeight); - // precompute interpolation data for box scheme for each cut element - if (isBox()) + // precompute interpolation data for box scheme for each cut bulk element + if constexpr (isBox()) { - const auto bulkElement = bulkGridGeometry.element(bulkElementIdx); - circleCornerIndices.push_back(&(this->vertexIndices(bulkIdx, bulkElementIdx))); - - // evaluate shape functions at the integration point - const auto bulkGeometry = bulkElement.geometry(); - ShapeValues shapeValues; - this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, circlePoints[k], shapeValues); - circleShapeValues.emplace_back(std::move(shapeValues)); + if (!static_cast(circleCornerIndices.count(bulkElementIdx))) + { + const auto bulkElement = this->problem(bulkIdx).fvGridGeometry().element(bulkElementIdx); + circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx); + + // evaluate shape functions at the integration point + const auto bulkGeometry = bulkElement.geometry(); + this->getShapeValues(bulkIdx, this->problem(bulkIdx).fvGridGeometry(), bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); + } } } } @@ -397,12 +398,10 @@ public: for (auto bulkElementIdx : bulkElementIndices) { const auto id = this->idCounter_++; - const auto ie = lowDimGeometry.integrationElement(qp.position()); - const auto qpweight = qp.weight(); - this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, ie, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({bulkElementIdx})); this->pointSources(bulkIdx).back().setEmbeddings(bulkElementIndices.size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({lowDimElementIdx})); this->pointSources(lowDimIdx).back().setEmbeddings(bulkElementIndices.size()); // pre compute additional data used for the evaluation of -- GitLab From 68f47d65b28b82f03798f5eecad0570395aba478 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 12 Feb 2019 18:10:37 +0100 Subject: [PATCH 04/17] [embedded] Use circle average for cylindersource mode --- .../embedded/couplingmanager1d3d.hh | 164 +++++++++++++++--- 1 file changed, 141 insertions(+), 23 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index fc4581ad85..fd3993147b 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -584,10 +584,11 @@ private: */ template class EmbeddedCouplingManager1d3d -: public EmbeddedCouplingManagerBase> +: public EmbeddedCouplingManagerBase, + CircleAveragePointSourceTraits> { using ThisType = EmbeddedCouplingManager1d3d; - using ParentType = EmbeddedCouplingManagerBase; + using ParentType = EmbeddedCouplingManagerBase>; using Scalar = typename MDTraits::Scalar; using SolutionVector = typename MDTraits::SolutionVector; using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData; @@ -627,6 +628,35 @@ public: computeLowDimVolumeFractions(); } + /*! + * \brief extend the jacobian pattern of the diagonal block of domain i + * by those entries that are not already in the uncoupled pattern + * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source + * term depends on a non-local average of a quantity of the same domain + */ + template + void extendJacobianPattern(Dune::index_constant domainI, JacobianPattern& pattern) const + { + extendedSourceStencil_.extendJacobianPattern(*this, domainI, pattern); + } + + /*! + * \brief evaluate additional derivatives of the element residual of a domain with respect + * to dofs in the same domain that are not in the regular stencil (per default this is not the case) + * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source + * term depends on a non-local average of a quantity of the same domain + * \note This is the same for box and cc + */ + template + void evalAdditionalDomainDerivatives(Dune::index_constant domainI, + const LocalAssemblerI& localAssemblerI, + const typename LocalAssemblerI::LocalResidual::ElementResidualVector&, + JacobianMatrixDiagBlock& A, + GridVariables& gridVariables) + { + extendedSourceStencil_.evalAdditionalDomainDerivatives(*this, domainI, localAssemblerI, this->curSol(), A, gridVariables); + } + /* \brief Compute integration point point sources and associated data * * This method uses grid glue to intersect the given grids. Over each intersection @@ -651,6 +681,7 @@ public: // clear all internal members like pointsource vectors and stencils // initializes the point source id counter this->clear(); + extendedSourceStencil_.stencil().clear(); // precompute the vertex indices for efficiency this->precomputeVertexIndices(bulkIdx); @@ -693,25 +724,81 @@ public: const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx); const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0); const auto integrationElement = lowDimGeometry.integrationElement(qp.position())*2*M_PI*radius/Scalar(numIp); - const auto weight = qp.weight()/(2*M_PI*radius); + const auto qpweight = qp.weight()/(2*M_PI*radius); + const auto circleAvgWeight = 2*M_PI*radius/numIp; const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); + std::vector> circleBulkElementIndices(circlePoints.size()); + std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); + std::vector circleStencil; circleStencil.reserve(circlePoints.size()); + // for box + std::unordered_map > circleCornerIndices; + using ShapeValues = std::vector >; + std::unordered_map circleShapeValues; + + // go over all circle points and precompute some quantities for the circle average + for (int k = 0; k < circlePoints.size(); ++k) + { + circleBulkElementIndices[k] = intersectingEntities(circlePoints[k], bulkTree); + if (circleBulkElementIndices[k].empty()) + continue; + + const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices[k].size(); + for (const auto bulkElementIdx : circleBulkElementIndices[k]) + { + circleStencil.push_back(bulkElementIdx); + circleIpWeight.push_back(localCircleAvgWeight); + + // precompute interpolation data for the box scheme + if (isBox()) + { + if (!static_cast(circleCornerIndices.count(bulkElementIdx))) + { + const auto bulkElement = this->problem(bulkIdx).fvGridGeometry().element(bulkElementIdx); + circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx); + + // evaluate shape functions at the integration point + const auto bulkGeometry = bulkElement.geometry(); + this->getShapeValues(bulkIdx, this->problem(bulkIdx).fvGridGeometry(), bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); + } + } + } + } + + // export low dim circle stencil + if (isBox()) + { + // we insert all vertices and make it unique later + for (const auto& vertices : circleCornerIndices) + { + this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), + vertices.second.begin(), vertices.second.end()); + + } + } + else + { + this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), + circleStencil.begin(), circleStencil.end()); + } + for (int k = 0; k < circlePoints.size(); ++k) { const auto& circlePos = circlePoints[k]; - const auto circleBulkElementIndices = intersectingEntities(circlePos, bulkTree); - if (circleBulkElementIndices.empty()) + + // find bulk elements intersection with the circle elements + if (circleBulkElementIndices[k].empty()) continue; // loop over the bulk elements at the integration points (usually one except when it is on a face or edge or vertex) // and add a point source at every point on the circle - for (const auto bulkElementIdx : circleBulkElementIndices) + for (const auto bulkElementIdx : circleBulkElementIndices[k]) { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(circlePos, id, weight, integrationElement, std::vector({bulkElementIdx})); - this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices.size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, weight, integrationElement, std::vector({lowDimElementIdx})); - this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices.size()); + this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); // pre compute additional data used for the evaluation of // the actual solution dependent source term @@ -732,6 +819,8 @@ public: // add data needed to compute integral over the circle if constexpr (isBox()) { + psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil); + using ShapeValues = std::vector >; const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); ShapeValues shapeValues; @@ -740,25 +829,13 @@ public: } else { + psData.addCircleInterpolation(circleIpWeight, circleStencil); psData.addBulkInterpolation(bulkElementIdx); } // publish point source data in the global vector this->pointSourceData().emplace_back(std::move(psData)); - // export the lowdim coupling stencil - // we insert all vertices / elements and make it unique later - if (isBox()) - { - const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); - this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), - vertices.begin(), vertices.end()); - } - else - { - this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); - } - // export the bulk coupling stencil // we insert all vertices / elements and make it unique later if (isBox()) @@ -773,11 +850,50 @@ public: this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); } + // export bulk circle stencil + if (isBox()) + { + // we insert all vertices and make it unique later + for (const auto& vertices : circleCornerIndices) + { + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + vertices.second.begin(), vertices.second.end()); + + } + } + else + { + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + circleStencil.begin(), circleStencil.end()); + } } } } } + // make the circle stencil unique (for source derivatives) + for (auto&& stencil : extendedSourceStencil_.stencil()) + { + std::sort(stencil.second.begin(), stencil.second.end()); + stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + + // remove the vertices element (box) + if (isBox()) + { + const auto& indices = this->vertexIndices(bulkIdx, stencil.first); + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), + stencil.second.end()); + } + // remove the own element (cell-centered) + else + { + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return i == stencil.first; }), + stencil.second.end()); + } + } + // make stencils unique using namespace Dune::Hybrid; forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) @@ -851,6 +967,8 @@ public: // \} private: + //! the extended source stencil object + EmbeddedCoupling::ExtendedSourceStencil extendedSourceStencil_; //! vector for the volume fraction of the lowdim domain in the bulk domain cells std::vector lowDimVolumeInBulkElement_; -- GitLab From d17cc232362d7c260b63dadf182a41f8af0013a9 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Mon, 4 Mar 2019 19:33:19 +0100 Subject: [PATCH 05/17] [temp] Improvements in the coupling manager --- .../embedded/couplingmanager1d3d.hh | 271 +++++++++--------- 1 file changed, 133 insertions(+), 138 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index fd3993147b..a784cfd225 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -34,6 +34,7 @@ #include #include +#include #include #include #include @@ -299,15 +300,20 @@ public: this->precomputeVertexIndices(bulkIdx); this->precomputeVertexIndices(lowDimIdx); - // iterate over all lowdim elements - const auto& lowDimProblem = this->problem(lowDimIdx); - for (const auto& lowDimElement : elements(this->gridView(lowDimIdx))) + // intersect the bounding box trees + this->glueGrids(); + + // iterate over all intersection and add point sources + for (const auto& is : intersections(this->glue())) { - // get the Gaussian quadrature rule for the low dim element - const auto lowDimGeometry = lowDimElement.geometry(); - const auto& quad = Dune::QuadratureRules::rule(lowDimGeometry.type(), order); + // all inside elements are identical... + const auto& lowDimElement = is.targetEntity(0); + const auto lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement); - const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(lowDimElement); + // get the intersection geometry + const auto intersectionGeometry = is.geometry(); + // get the Gaussian quadrature rule for the local intersection + const auto& quad = Dune::QuadratureRules::rule(intersectionGeometry.type(), order); // apply the Gaussian quadrature rule and define point sources at each quadrature point // note that the approximation is not optimal if @@ -319,7 +325,7 @@ public: for (auto&& qp : quad) { // global position of the quadrature point - const auto globalPos = lowDimGeometry.global(qp.position()); + const auto globalPos = intersectionGeometry.global(qp.position()); const auto bulkElementIndices = intersectingEntities(globalPos, bulkTree); @@ -333,10 +339,10 @@ public: ////////////////////////////////////////////////////////// static const auto numIp = getParam("MixedDimension.NumCircleSegments"); - const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx); - const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0); + const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); + const auto normal = intersectionGeometry.corner(1)-intersectionGeometry.corner(0); const auto circleAvgWeight = 2*M_PI*radius/numIp; - const auto integrationElement = lowDimGeometry.integrationElement(qp.position()); + const auto integrationElement = intersectionGeometry.integrationElement(qp.position()); const auto qpweight = qp.weight(); const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); @@ -366,12 +372,12 @@ public: { if (!static_cast(circleCornerIndices.count(bulkElementIdx))) { - const auto bulkElement = this->problem(bulkIdx).fvGridGeometry().element(bulkElementIdx); + const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx); // evaluate shape functions at the integration point const auto bulkGeometry = bulkElement.geometry(); - this->getShapeValues(bulkIdx, this->problem(bulkIdx).fvGridGeometry(), bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); } } } @@ -411,7 +417,7 @@ public: if constexpr (isBox()) { ShapeValues shapeValues; - this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues); + this->getShapeValues(lowDimIdx, lowDimFvGridGeometry, intersectionGeometry, globalPos, shapeValues); psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); } else @@ -424,9 +430,9 @@ public: { psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil); - const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); + const auto bulkGeometry = bulkFvGridGeometry.element(bulkElementIdx).geometry(); ShapeValues shapeValues; - this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, globalPos, shapeValues); + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, globalPos, shapeValues); psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); } else @@ -438,6 +444,21 @@ public: // publish point source data in the global vector this->pointSourceData().emplace_back(std::move(psData)); + // compute the average flux scaling factor + const auto outsideGeometry = bulkFvGridGeometry.element(bulkElementIdx).geometry(); + const auto& quadBulk = Dune::QuadratureRules::rule(outsideGeometry.type(), 5); + // iterate over all quadrature points + Scalar avgMinDist = 0.0; + for (auto&& qpBulk : quadBulk) + { + const auto globalPosBulk = outsideGeometry.global(qpBulk.position()); + const auto d = computeMinDistance_(intersectionGeometry.corner(0), intersectionGeometry.corner(1), globalPosBulk); + avgMinDist += d*qpBulk.weight(); + } + + // add it to the internal data vector + this->averageDistanceToBulkCell().push_back(avgMinDist); + // export the bulk coupling stencil if (isBox()) { @@ -567,6 +588,16 @@ public: // \} private: + inline Scalar computeMinDistance_(const GlobalPosition& a, const GlobalPosition& b, const GlobalPosition& x3d) const + { + // r is the minimum distance to the line through a and b + const auto ab = b - a; + const auto t = (x3d - a)*ab/ab.two_norm2(); + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - x3d).two_norm(); + return r; + } + //! the extended source stencil object EmbeddedCoupling::ExtendedSourceStencil extendedSourceStencil_; @@ -728,13 +759,13 @@ public: const auto circleAvgWeight = 2*M_PI*radius/numIp; const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); - std::vector> circleBulkElementIndices(circlePoints.size()); + std::vector>> circleBulkElementIndices(circlePoints.size()); std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); - std::vector circleStencil; circleStencil.reserve(circlePoints.size()); + std::vector> circleStencil; circleStencil.reserve(circlePoints.size()); // for box - std::unordered_map > circleCornerIndices; + std::unordered_map, std::vector> > circleCornerIndices; using ShapeValues = std::vector >; - std::unordered_map circleShapeValues; + std::unordered_map, ShapeValues> circleShapeValues; // go over all circle points and precompute some quantities for the circle average for (int k = 0; k < circlePoints.size(); ++k) @@ -1097,13 +1128,26 @@ public: this->averageDistanceToBulkCell().reserve(this->glue().size()); fluxScalingFactor_.reserve(this->glue().size()); + // reserve memory for stencils + const auto numBulkElements = this->gridView(bulkIdx).size(0); + for (GridIndex bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) + { + this->couplingStencils(bulkIdx)[bulkElementIdx].reserve(50); + extendedSourceStencil_.stencil()[bulkElementIdx].reserve(50); + bulkSourceIds_[bulkElementIdx][0].reserve(20); + bulkSourceWeights_[bulkElementIdx][0].reserve(20); + } + // generate a bunch of random vectors and values for // Monte-carlo integration on the cylinder defined by line and radius - static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); - const int samples = useCylinderIntegration ? getParam("MixedDimension.NumMCSamples") : 1; - CylinderIntegration cylIntegration(samples); + static const int samples = getParam("MixedDimension.NumMCSamples"); + CylinderIntegration cylIntegration(samples, 1); static const auto kernelWidthFactor = getParam("MixedDimension.KernelWidthFactor"); + std::ofstream file("points.log", std::ios::app); + file << "x, y, z\n"; + file.close(); + for (const auto& is : intersections(this->glue())) { // all inside elements are identical... @@ -1118,12 +1162,12 @@ public: // * 3d: a new kernel volume source const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); const auto kernelWidth = kernelWidthFactor*radius; - const auto a = intersectionGeometry.corner(0); const auto b = intersectionGeometry.corner(1); + cylIntegration.setGeometry(a, b, kernelWidth); // we can have multiple 3d elements per 1d intersection, for evaluation taking one of them should be fine - for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) + for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) { // compute source id // for each id we have @@ -1139,11 +1183,9 @@ public: this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); const auto& outside = is.domainEntity(outsideIdx); - const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside); + const auto bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside); // compute the weights for the bulk volume sources - if (useCylinderIntegration) - cylIntegration.setGeometry(a, b, kernelWidth); computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); // pre compute additional data used for the evaluation of @@ -1310,11 +1352,11 @@ public: } //! return all source ids for a bulk elements - const std::vector& bulkSourceIds(GridIndex eIdx, std::size_t scvIdx = 0) const + const std::vector& bulkSourceIds(GridIndex eIdx, int scvIdx = 0) const { return bulkSourceIds_[eIdx][scvIdx]; } //! return all source ids for a bulk elements - const std::vector& bulkSourceWeights(GridIndex, std::size_t scvIdx = 0) const + const std::vector& bulkSourceWeights(GridIndex eIdx, int scvIdx = 0) const { return bulkSourceWeights_[eIdx][scvIdx]; } //! The flux scaling factor for a source with id @@ -1327,111 +1369,45 @@ private: //! compute the kernel distributed sources and add stencils template void computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, - std::size_t id, std::size_t lowDimElementIdx, std::size_t coupledBulkElementIdx, - const CylIntegration& cylIntegration, std::size_t embeddings) + std::size_t id, GridIndex lowDimElementIdx, GridIndex coupledBulkElementIdx, + const CylIntegration& cylIntegration, int embeddings) { // Monte-carlo integration on the cylinder defined by line and radius - const auto numBulkElements = this->gridView(bulkIdx).size(0); - constexpr std::size_t numScv = isBox() ? 1<> weights(numBulkElements, std::array{}); - std::vector> visited(numBulkElements, std::bitset{}); + static const auto min = getParam("Tissue.Grid.LowerLeft"); + static const auto max = getParam("Tissue.Grid.UpperRight"); + static const auto cells = getParam("Tissue.Grid.Cells"); + const auto cylSamples = cylIntegration.size(); + const auto& a = line.corner(0); + const auto& b = line.corner(1); - Scalar integral = 0.0; - static const bool useCylinderIntegration = getParam("MixedDimension.UseCylinderIntegration"); - if (useCylinderIntegration) - { - for (int i = 0; i < cylIntegration.size(); ++i) - { - const auto point = cylIntegration.getIntegrationPoint(i); - const auto weight = evalKernel(line.corner(0), line.corner(1), point, radius, kernelWidth)*cylIntegration.integrationElement(i); - integral += weight; - - static const auto min = getParam("Tissue.Grid.LowerLeft"); - static const auto max = getParam("Tissue.Grid.UpperRight"); - static const auto cells = getParam("Tissue.Grid.Cells"); - const auto is = intersectingEntityCartesian(point, min, max, cells); - // const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); - if (is.first) - { - const auto bulkElementIdx = is.second; - if (isBox()) - { - const auto dist = point-min; - const auto diag = max-min; - std::bitset<3> corner; - for (int i = 0; i < 3; ++i) - corner[i] = bool(int(std::round(dist[i]/diag[i]))); - - std::size_t scvIdx = corner.to_ulong(); - visited[bulkElementIdx][scvIdx] = true; - weights[bulkElementIdx][scvIdx] += weight; - } - else - { - visited[bulkElementIdx][0] = true; - weights[bulkElementIdx][0] += weight; - } - } - } - } - // integrate over 3d domain using standard quadrature rules - else + std::ofstream file("points.log", std::ios::app); + + for (int i = 0; i < cylSamples; ++i) { - static const int maxLevel = getParam("MixedDimension.MaxOctreeLevel"); - CylinderDomain cylinder(line.corner(0), line.corner(1), kernelWidth); - const auto kernelFunc = [&](const auto& p){ return evalKernel(line.corner(0), line.corner(1), p, radius, kernelWidth); }; - for (const auto& element : elements(this->gridView(bulkIdx))) + const auto& point = cylIntegration.getIntegrationPoint(i); + file << point[0] << ", " << point[1] << ", " << point[2] << "\n"; + //const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); + if (const auto [hasIntersection, bulkElementIdx] = intersectingEntityCartesian(point, min, max, cells); hasIntersection) { - const auto bulkElementIdx = this->problem(bulkIdx).fvGridGeometry().elementMapper().index(element); - const auto geometry = element.geometry(); - OctreeIntegrator integrator(geometry.corner(0), geometry.corner(7), maxLevel); - const auto weight = integrator.integrate(kernelFunc, cylinder); - if (weight > 1e-35) + const auto localWeight = evalConstKernel_(a, b, point, radius, kernelWidth)*cylIntegration.integrationElement(i)/embeddings; + if (!bulkSourceIds_[bulkElementIdx][0].empty() && id == bulkSourceIds_[bulkElementIdx][0].back()) { - visited[bulkElementIdx][0] = true; - weights[bulkElementIdx][0] += weight; - integral += weight; + bulkSourceWeights_[bulkElementIdx][0].back() += localWeight; } - } - } - - // const auto length = line.volume(); - // std::cout << "integration error in %: " << std::abs(length-integral)/length*100 << std::endl; - // const auto correctionFactor = length/integral; - for (const auto& element : elements(this->problem(bulkIdx).fvGridGeometry().gridView())) - { - const auto bulkElementIdx = this->problem(bulkIdx).fvGridGeometry().elementMapper().index(element); - if (!isBox() && !visited[bulkElementIdx][0]) - continue; - - if (isBox()) - { - auto fvGeometry = localView(this->problem(bulkIdx).fvGridGeometry()); - fvGeometry.bindElement(element); - for (const auto& scv : scvs(fvGeometry)) + else { - const auto scvIdx = scv.indexInElement(); - if (!visited[bulkElementIdx][scvIdx]) - continue; - - bulkSourceIds_[bulkElementIdx][scvIdx].push_back(id); - bulkSourceWeights_[bulkElementIdx][scvIdx].push_back(weights[bulkElementIdx][scvIdx]/embeddings);//*correctionFactor); + bulkSourceIds_[bulkElementIdx][0].emplace_back(id); + bulkSourceWeights_[bulkElementIdx][0].emplace_back(localWeight); addBulkSourceStencils_(bulkElementIdx, lowDimElementIdx, coupledBulkElementIdx); } - - } - else - { - bulkSourceIds_[bulkElementIdx][0].push_back(id); - bulkSourceWeights_[bulkElementIdx][0].push_back(weights[bulkElementIdx][0]/embeddings);//*correctionFactor); - addBulkSourceStencils_(bulkElementIdx, lowDimElementIdx, coupledBulkElementIdx); } } + } //! add additional stencil entries for the bulk element - void addBulkSourceStencils_(std::size_t bulkElementIdx, std::size_t coupledLowDimElementIdx, std::size_t coupledBulkElementIdx) + void addBulkSourceStencils_(GridIndex bulkElementIdx, GridIndex coupledLowDimElementIdx, GridIndex coupledBulkElementIdx) { // add the lowdim element to the coupling stencil of this bulk element if (isBox()) @@ -1443,7 +1419,8 @@ private: } else { - this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(coupledLowDimElementIdx); + auto& s = this->couplingStencils(bulkIdx)[bulkElementIdx]; + s.push_back(coupledLowDimElementIdx); } // the extended source stencil, every 3d element with a source is coupled to @@ -1456,16 +1433,17 @@ private: } else { - extendedSourceStencil_.stencil()[bulkElementIdx].push_back(coupledBulkElementIdx); + auto& s = extendedSourceStencil_.stencil()[bulkElementIdx]; + s.push_back(coupledBulkElementIdx); } } //! a cylindrical kernel around the segment a->b - inline Scalar evalKernel(const GlobalPosition& a, - const GlobalPosition& b, - const GlobalPosition& point, - const Scalar R, - const Scalar rho) const + Scalar evalConstKernel_(const GlobalPosition& a, + const GlobalPosition& b, + const GlobalPosition& point, + const Scalar R, + const Scalar rho) const noexcept { // projection of point onto line a + t*(b-a) const auto ab = b - a; @@ -1482,20 +1460,37 @@ private: if (r > rho) return 0.0; - static const std::string type = getParam("MixedDimension.KernelType"); - if (type == "Const") - return 1.0/(M_PI*rho*rho); - else if (type == "Cubic") - { - const auto r2 = r*r; const auto r3 = r*r*r; - const auto rho2 = rho*rho; const auto rho3 = rho2*rho; - return 10.0*(2.0*r3/rho3 - 3.0*r2/rho2 + 1.0)/(3.0*M_PI*rho2); - } + return 1.0/(M_PI*rho*rho); + } + + //! a cylindrical kernel around the segment a->b + Scalar evalCubicKernel_(const GlobalPosition& a, + const GlobalPosition& b, + const GlobalPosition& point, + const Scalar R, + const Scalar rho) const noexcept + { + // projection of point onto line a + t*(b-a) + const auto ab = b - a; + const auto t = (point - a)*ab/ab.two_norm2(); - DUNE_THROW(Dune::NotImplemented, "Kernel type: " << type); + // return 0 if we are outside cylinder + if (t < 0.0 || t > 1.0) + return 0.0; + + // compute distance + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - point).two_norm(); + + if (r > rho) + return 0.0; + + const auto r2 = r*r; const auto r3 = r*r*r; + const auto rho2 = rho*rho; const auto rho3 = rho2*rho; + return 10.0*(2.0*r3/rho3 - 3.0*r2/rho2 + 1.0)/(3.0*M_PI*rho2); } - inline Scalar computeMinDistance_(const GlobalPosition& a, const GlobalPosition& b, const GlobalPosition& x3d) const + Scalar computeMinDistance_(const GlobalPosition& a, const GlobalPosition& b, const GlobalPosition& x3d) const noexcept { // r is the minimum distance to the line through a and b const auto ab = b - a; @@ -1510,7 +1505,7 @@ private: //! vector for the volume fraction of the lowdim domain in the bulk domain cells std::vector lowDimVolumeInBulkElement_; //! kernel sources to integrate for each bulk element - std::vector, isBox() ? 1<> bulkSourceIds_; + std::vector, isBox() ? 1<> bulkSourceIds_; //! the integral of the kernel for each point source / integration point, i.e. weight for the source std::vector, isBox() ? 1<> bulkSourceWeights_; //! the flux scaling factor for the respective source id -- GitLab From 4d507597a6cf7b490eacb667c700c3262dafac17 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Tue, 5 Mar 2019 11:57:40 +0100 Subject: [PATCH 06/17] [1d3d] Surface factor --- .../embedded/couplingmanager1d3d.hh | 49 ++++++++++++------- 1 file changed, 30 insertions(+), 19 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index a784cfd225..ba8d809c81 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -361,6 +361,7 @@ public: if (circleBulkElementIndices.empty()) continue; + ++insideCirclePoints; const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices.size(); for (const auto bulkElementIdx : circleBulkElementIndices) { @@ -400,14 +401,16 @@ public: circleStencil.begin(), circleStencil.end()); } + // surface fraction that is inside the domain (only different than 1.0 on the boundary) + const auto surfaceFraction = Scalar(insideCirclePoints)/Scalar(circlePoints.size()); // loop over the bulk elements at the integration points (usually one except when it is on a face or edge or vertex) for (auto bulkElementIdx : bulkElementIndices) { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, std::vector({bulkElementIdx})); this->pointSources(bulkIdx).back().setEmbeddings(bulkElementIndices.size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, std::vector({lowDimElementIdx})); this->pointSources(lowDimIdx).back().setEmbeddings(bulkElementIndices.size()); // pre compute additional data used for the evaluation of @@ -768,12 +771,14 @@ public: std::unordered_map, ShapeValues> circleShapeValues; // go over all circle points and precompute some quantities for the circle average + int insideCirclePoints = 0; for (int k = 0; k < circlePoints.size(); ++k) { circleBulkElementIndices[k] = intersectingEntities(circlePoints[k], bulkTree); if (circleBulkElementIndices[k].empty()) continue; + ++insideCirclePoints; const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices[k].size(); for (const auto bulkElementIdx : circleBulkElementIndices[k]) { @@ -813,6 +818,9 @@ public: circleStencil.begin(), circleStencil.end()); } + // surface fraction that is inside the domain (only different than 1.0 on the boundary) + const auto surfaceFraction = Scalar(insideCirclePoints)/Scalar(circlePoints.size()); + for (int k = 0; k < circlePoints.size(); ++k) { const auto& circlePos = circlePoints[k]; @@ -826,9 +834,10 @@ public: for (const auto bulkElementIdx : circleBulkElementIndices[k]) { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement, std::vector({bulkElementIdx})); + + this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement*surfaceFraction, std::vector({bulkElementIdx})); this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, std::vector({lowDimElementIdx})); this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); // pre compute additional data used for the evaluation of @@ -1177,16 +1186,15 @@ public: // (3) the flux scaling factor for each outside element, i.e. each id const auto id = this->idCounter_++; - // place a point source at the intersection center - const auto center = intersectionGeometry.center(); - this->pointSources(lowDimIdx).emplace_back(center, id, /*weight=*/1.0, intersectionGeometry.volume(), std::vector({lowDimElementIdx})); - this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); - + // compute the weights for the bulk volume sources const auto& outside = is.domainEntity(outsideIdx); const auto bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside); + const auto surfaceFactor = computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); - // compute the weights for the bulk volume sources - computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); + // place a point source at the intersection center + const auto center = intersectionGeometry.center(); + this->pointSources(lowDimIdx).emplace_back(center, id, surfaceFactor, intersectionGeometry.volume(), std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); // pre compute additional data used for the evaluation of // the actual solution dependent source term @@ -1368,9 +1376,9 @@ public: private: //! compute the kernel distributed sources and add stencils template - void computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, - std::size_t id, GridIndex lowDimElementIdx, GridIndex coupledBulkElementIdx, - const CylIntegration& cylIntegration, int embeddings) + Scalar computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, + std::size_t id, GridIndex lowDimElementIdx, GridIndex coupledBulkElementIdx, + const CylIntegration& cylIntegration, int embeddings) { // Monte-carlo integration on the cylinder defined by line and radius static const auto min = getParam("Tissue.Grid.LowerLeft"); @@ -1380,16 +1388,17 @@ private: const auto& a = line.corner(0); const auto& b = line.corner(1); - std::ofstream file("points.log", std::ios::app); - + //std::ofstream file("points.log", std::ios::app); + Scalar integral = 0.0; for (int i = 0; i < cylSamples; ++i) { const auto& point = cylIntegration.getIntegrationPoint(i); - file << point[0] << ", " << point[1] << ", " << point[2] << "\n"; + //file << point[0] << ", " << point[1] << ", " << point[2] << "\n"; //const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); if (const auto [hasIntersection, bulkElementIdx] = intersectingEntityCartesian(point, min, max, cells); hasIntersection) { - const auto localWeight = evalConstKernel_(a, b, point, radius, kernelWidth)*cylIntegration.integrationElement(i)/embeddings; + const auto localWeight = evalConstKernel_(a, b, point, radius, kernelWidth)*cylIntegration.integrationElement(i)/Scalar(embeddings); + integral += localWeight; if (!bulkSourceIds_[bulkElementIdx][0].empty() && id == bulkSourceIds_[bulkElementIdx][0].back()) { bulkSourceWeights_[bulkElementIdx][0].back() += localWeight; @@ -1403,7 +1412,9 @@ private: } } - + // the surface factor (which fraction of the source is inside the domain and needs to be considered) + const auto length = (a-b).two_norm()/Scalar(embeddings); + return integral/length; } //! add additional stencil entries for the bulk element -- GitLab From 140aae25a919abc6c637169aa8b3edbfe77bccf6 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 15 Mar 2019 20:32:51 +0100 Subject: [PATCH 07/17] [1d3d] Add to string function for embedded coupling modes --- .../multidomain/embedded/couplingmanager1d3d.hh | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index ba8d809c81..4a436e0e78 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -53,6 +53,22 @@ namespace Dumux { enum class EmbeddedCouplingMode { line, average, cylindersources, kernel }; +/*! + * \ingroup EmbeddedCoupling + * \brief Convert the coupling mode to a human readable name + */ +std::string toString(EmbeddedCouplingMode mode) +{ + switch (mode) + { + case EmbeddedCouplingMode::line: return "line"; + case EmbeddedCouplingMode::average: return "average"; + case EmbeddedCouplingMode::cylindersources: return "cylinder"; + case EmbeddedCouplingMode::kernel: return "kernel"; + default: return "unknown"; + } +} + //! point source traits for the circle average coupling mode template struct CircleAveragePointSourceTraits -- GitLab From 275d3caa5d8b7d2c38720c65ae402710d82bc172 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 15 Mar 2019 20:33:34 +0100 Subject: [PATCH 08/17] [1d3d] Fix cylinder coupling weights --- .../embedded/couplingmanager1d3d.hh | 33 +++++++++++-------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 4a436e0e78..631a35de27 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -737,15 +737,22 @@ public: this->precomputeVertexIndices(bulkIdx); this->precomputeVertexIndices(lowDimIdx); + // intersect the bounding box trees + this->glueGrids(); + // iterate over all lowdim elements const auto& lowDimProblem = this->problem(lowDimIdx); - for (const auto& lowDimElement : elements(this->gridView(lowDimIdx))) + // iterate over all intersection and add point sources + for (const auto& is : intersections(this->glue())) { - // get the Gaussian quadrature rule for the low dim element - const auto lowDimGeometry = lowDimElement.geometry(); - const auto& quad = Dune::QuadratureRules::rule(lowDimGeometry.type(), order); + // all inside elements are identical... + const auto& lowDimElement = is.inside(0); + const auto lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement); - const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(lowDimElement); + // get the intersection geometry + const auto intersectionGeometry = is.geometry(); + // get the Gaussian quadrature rule for the local intersection + const auto& quad = Dune::QuadratureRules::rule(intersectionGeometry.type(), order); // apply the Gaussian quadrature rule and define point sources at each quadrature point // note that the approximation is not optimal if @@ -757,7 +764,7 @@ public: for (auto&& qp : quad) { // global position of the quadrature point - const auto globalPos = lowDimGeometry.global(qp.position()); + const auto globalPos = intersectionGeometry.global(qp.position()); const auto bulkElementIndices = intersectingEntities(globalPos, bulkTree); @@ -770,15 +777,15 @@ public: // get points on the cylinder surface at the integration point //////////////////////////////////////////////////////////////// - static const auto numIp = getParam("MixedDimension.NumCircleSegments", 25); + static const auto numIp = getParam("MixedDimension.NumCircleSegments"); const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx); - const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0); - const auto integrationElement = lowDimGeometry.integrationElement(qp.position())*2*M_PI*radius/Scalar(numIp); + const auto normal = intersectionGeometry.corner(1)-intersectionGeometry.corner(0); + const auto integrationElement = intersectionGeometry.integrationElement(qp.position())*2*M_PI*radius/Scalar(numIp); const auto qpweight = qp.weight()/(2*M_PI*radius); const auto circleAvgWeight = 2*M_PI*radius/numIp; const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); - std::vector>> circleBulkElementIndices(circlePoints.size()); + std::vector> circleBulkElementIndices(circlePoints.size()); std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); std::vector> circleStencil; circleStencil.reserve(circlePoints.size()); // for box @@ -806,12 +813,12 @@ public: { if (!static_cast(circleCornerIndices.count(bulkElementIdx))) { - const auto bulkElement = this->problem(bulkIdx).fvGridGeometry().element(bulkElementIdx); + const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx); // evaluate shape functions at the integration point const auto bulkGeometry = bulkElement.geometry(); - this->getShapeValues(bulkIdx, this->problem(bulkIdx).fvGridGeometry(), bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); } } } @@ -864,7 +871,7 @@ public: { using ShapeValues = std::vector >; ShapeValues shapeValues; - this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues); + this->getShapeValues(lowDimIdx, lowDimGridGeometry, intersectionGeometry, globalPos, shapeValues); psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); } else -- GitLab From 23b214bb9c6e3313a74de0995c3f94436b9f2c68 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 22 Mar 2019 13:32:34 +0100 Subject: [PATCH 09/17] [1d3d] Add projection method to coupling modes --- dumux/multidomain/embedded/couplingmanager1d3d.hh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 631a35de27..404fc0951d 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -51,7 +51,7 @@ namespace Dumux { * \brief The coupling mode */ enum class EmbeddedCouplingMode -{ line, average, cylindersources, kernel }; +{ line, average, cylindersources, kernel, projection }; /*! * \ingroup EmbeddedCoupling @@ -65,6 +65,7 @@ std::string toString(EmbeddedCouplingMode mode) case EmbeddedCouplingMode::average: return "average"; case EmbeddedCouplingMode::cylindersources: return "cylinder"; case EmbeddedCouplingMode::kernel: return "kernel"; + case EmbeddedCouplingMode::projection: return "projection"; default: return "unknown"; } } -- GitLab From ac5666c6f9fa55379e92077e1c3cf968097b42b8 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 22 Mar 2019 13:33:41 +0100 Subject: [PATCH 10/17] [1d3d][fixup] Remove wrongly introduces surface factor for cylinder method --- dumux/multidomain/embedded/couplingmanager1d3d.hh | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 404fc0951d..c5c54477bf 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -795,14 +795,12 @@ public: std::unordered_map, ShapeValues> circleShapeValues; // go over all circle points and precompute some quantities for the circle average - int insideCirclePoints = 0; for (int k = 0; k < circlePoints.size(); ++k) { circleBulkElementIndices[k] = intersectingEntities(circlePoints[k], bulkTree); if (circleBulkElementIndices[k].empty()) continue; - ++insideCirclePoints; const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices[k].size(); for (const auto bulkElementIdx : circleBulkElementIndices[k]) { @@ -842,9 +840,6 @@ public: circleStencil.begin(), circleStencil.end()); } - // surface fraction that is inside the domain (only different than 1.0 on the boundary) - const auto surfaceFraction = Scalar(insideCirclePoints)/Scalar(circlePoints.size()); - for (int k = 0; k < circlePoints.size(); ++k) { const auto& circlePos = circlePoints[k]; @@ -859,9 +854,9 @@ public: { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement*surfaceFraction, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement, std::vector({bulkElementIdx})); this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({lowDimElementIdx})); this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); // pre compute additional data used for the evaluation of -- GitLab From dabfed82ff90023b18e61a321109d1137336f1ab Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 22 Mar 2019 13:34:36 +0100 Subject: [PATCH 11/17] [1d3d] Add restriction for box method --- dumux/multidomain/embedded/couplingmanager1d3d.hh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index c5c54477bf..dab169c1c2 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -1068,6 +1068,8 @@ class EmbeddedCouplingManager1d3d static constexpr bool isBox() { return GridGeometry::discMethod == DiscretizationMethod::box; } + static_assert(!isBox() && !isBox(), "The kernel coupling method is only implemented for the tpfa method"); + enum { bulkDim = GridView::dimension, lowDimDim = GridView::dimension, -- GitLab From c084d60ec450e8b5c665b34746d8d58f2a27cd2b Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 29 Mar 2019 15:40:38 +0100 Subject: [PATCH 12/17] [1d3d] Allow variable kernel width --- .../embedded/couplingmanager1d3d.hh | 60 ++++++++++++++++--- 1 file changed, 53 insertions(+), 7 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index dab169c1c2..78f3cba95e 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -1075,6 +1075,15 @@ class EmbeddedCouplingManager1d3d lowDimDim = GridView::dimension, dimWorld = GridView::dimensionworld }; + + // detect if a class (the spatial params) has a kernelWidthFactor() function + template + using VariableKernelWidthDetector = decltype(std::declval().kernelWidthFactor(std::declval()...)); + + template + static constexpr bool hasKernelWidthFactor() + { return Dune::Std::is_detected::value; } + public: static constexpr EmbeddedCouplingMode couplingMode = EmbeddedCouplingMode::kernel; @@ -1172,11 +1181,13 @@ public: // Monte-carlo integration on the cylinder defined by line and radius static const int samples = getParam("MixedDimension.NumMCSamples"); CylinderIntegration cylIntegration(samples, 1); - static const auto kernelWidthFactor = getParam("MixedDimension.KernelWidthFactor"); - std::ofstream file("points.log", std::ios::app); - file << "x, y, z\n"; - file.close(); + static const auto writeIntegrationPointsToFile = getParam("MixedDimension.WriteIntegrationPointsToFile", false); + if (writeIntegrationPointsToFile) + { + std::ofstream ipPointFile("kernel_points.log", std::ios::trunc); + ipPointFile << "x,y,z\n"; + } for (const auto& is : intersections(this->glue())) { @@ -1191,6 +1202,7 @@ public: // * 1d: a new point source // * 3d: a new kernel volume source const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); + const auto kernelWidthFactor = kernelWidthFactor_(this->problem(lowDimIdx).spatialParams(), lowDimElementIdx); const auto kernelWidth = kernelWidthFactor*radius; const auto a = intersectionGeometry.corner(0); const auto b = intersectionGeometry.corner(1); @@ -1267,7 +1279,7 @@ public: // add it to the internal data vector this->averageDistanceToBulkCell().push_back(avgMinDist); - fluxScalingFactor_.push_back(this->problem(bulkIdx).fluxScalingFactor(avgMinDist, radius)); + fluxScalingFactor_.push_back(this->problem(bulkIdx).fluxScalingFactor(avgMinDist, radius, kernelWidth)); // export the lowdim coupling stencil // we insert all vertices / elements and make it unique later @@ -1409,15 +1421,30 @@ private: const auto& a = line.corner(0); const auto& b = line.corner(1); - //std::ofstream file("points.log", std::ios::app); + // optionally write debugging / visual output of the integration points + static const auto writeIntegrationPointsToFile = getParam("MixedDimension.WriteIntegrationPointsToFile", false); + if (writeIntegrationPointsToFile) + { + std::ofstream ipPointFile("kernel_points.log", std::ios::app); + for (int i = 0; i < cylSamples; ++i) + { + const auto& point = cylIntegration.getIntegrationPoint(i); + const auto [hasIntersection, bulkElementIdx] = intersectingEntityCartesian(point, min, max, cells); + if (hasIntersection) + ipPointFile << point[0] << "," << point[1] << "," << point[2] << '\n'; + } + } + Scalar integral = 0.0; for (int i = 0; i < cylSamples; ++i) { const auto& point = cylIntegration.getIntegrationPoint(i); - //file << point[0] << ", " << point[1] << ", " << point[2] << "\n"; + // TODO: below only works for Cartesian grids with ijk numbering (e.g. level 0 YaspGrid (already fails when refined)) + // more general is the bounding box tree solution which always works, however it's much slower //const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); if (const auto [hasIntersection, bulkElementIdx] = intersectingEntityCartesian(point, min, max, cells); hasIntersection) { + assert(bulkElementIdx < this->problem(bulkIdx).fvGridGeometry().gridView().size(0)); const auto localWeight = evalConstKernel_(a, b, point, radius, kernelWidth)*cylIntegration.integrationElement(i)/Scalar(embeddings); integral += localWeight; if (!bulkSourceIds_[bulkElementIdx][0].empty() && id == bulkSourceIds_[bulkElementIdx][0].back()) @@ -1532,6 +1559,25 @@ private: return r; } + /*! + * \brief Get the kernel width factor from the spatial params (if possible) + */ + template + auto kernelWidthFactor_(const SpatialParams& spatialParams, unsigned int eIdx) + -> std::enable_if_t(), Scalar> + { return spatialParams.kernelWidthFactor(eIdx); } + + /*! + * \brief Get the kernel width factor (constant) from the input file (if possible) + */ + template + auto kernelWidthFactor_(const SpatialParams& spatialParams, unsigned int eIdx) + -> std::enable_if_t(), Scalar> + { + static const Scalar kernelWidthFactor = getParam("MixedDimension.KernelWidthFactor"); + return kernelWidthFactor; + } + //! the extended source stencil object for kernel coupling EmbeddedCoupling::ExtendedSourceStencil extendedSourceStencil_; //! vector for the volume fraction of the lowdim domain in the bulk domain cells -- GitLab From d7d94b85b928084827e09fd556977b04dd4d85fb Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 17 May 2019 19:53:26 +0200 Subject: [PATCH 13/17] [1d3d] Add static asserts for kernel method --- dumux/multidomain/embedded/couplingmanager1d3d.hh | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 78f3cba95e..35741b3cc1 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -1069,6 +1069,7 @@ class EmbeddedCouplingManager1d3d { return GridGeometry::discMethod == DiscretizationMethod::box; } static_assert(!isBox() && !isBox(), "The kernel coupling method is only implemented for the tpfa method"); + static_assert(Dune::Capabilities::isCartesian::Grid>::v, "The kernel coupling method is only implemented for structured grids"); enum { bulkDim = GridView::dimension, @@ -1095,6 +1096,10 @@ public: { ParentType::init(bulkProblem, lowDimProblem, curSol); computeLowDimVolumeFractions(); + + const auto refinement = getParamFromGroup(bulkProblem->paramGroup(), "Grid.Refinement", 0); + if (refinement > 0) + DUNE_THROW(Dune::NotImplemented, "The current intersection detection may likely fail for refined grids."); } /*! @@ -1414,9 +1419,10 @@ private: const CylIntegration& cylIntegration, int embeddings) { // Monte-carlo integration on the cylinder defined by line and radius - static const auto min = getParam("Tissue.Grid.LowerLeft"); - static const auto max = getParam("Tissue.Grid.UpperRight"); - static const auto cells = getParam("Tissue.Grid.Cells"); + static const auto bulkParamGroup = this->problem(bulkIdx).paramGroup(); + static const auto min = getParamFromGroup(bulkParamGroup, "Grid.LowerLeft"); + static const auto max = getParamFromGroup(bulkParamGroup, "Grid.UpperRight"); + static const auto cells = getParamFromGroup(bulkParamGroup, "Grid.Cells"); const auto cylSamples = cylIntegration.size(); const auto& a = line.corner(0); const auto& b = line.corner(1); -- GitLab From 31234c3bc8032fc59e3bee5da9c4f3540f6bf0be Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 17 May 2019 19:54:28 +0200 Subject: [PATCH 14/17] [1d3d] Remove deprecation warnings after point source update --- dumux/multidomain/embedded/couplingmanager1d3d.hh | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 35741b3cc1..cc81205c96 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -425,9 +425,9 @@ public: { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, bulkElementIdx); this->pointSources(bulkIdx).back().setEmbeddings(bulkElementIndices.size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, lowDimElementIdx); this->pointSources(lowDimIdx).back().setEmbeddings(bulkElementIndices.size()); // pre compute additional data used for the evaluation of @@ -854,9 +854,9 @@ public: { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement, bulkElementIdx); this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, lowDimElementIdx); this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); // pre compute additional data used for the evaluation of @@ -1231,7 +1231,7 @@ public: // place a point source at the intersection center const auto center = intersectionGeometry.center(); - this->pointSources(lowDimIdx).emplace_back(center, id, surfaceFactor, intersectionGeometry.volume(), std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(center, id, surfaceFactor, intersectionGeometry.volume(), lowDimElementIdx); this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); // pre compute additional data used for the evaluation of -- GitLab From 7f7180a8d8e32788d16e6b42257246827f8de2c4 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 10 Jul 2019 15:39:44 +0200 Subject: [PATCH 15/17] [1d3d] Slightly restructure kernel coupling manager internally --- .../embedded/couplingmanager1d3d.hh | 141 ++++++++++-------- 1 file changed, 77 insertions(+), 64 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index cc81205c96..a477bdba3a 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -1147,41 +1147,12 @@ public: Dune::Timer watch; std::cout << "Initializing the point sources..." << std::endl; - // clear all internal members like pointsource vectors and stencils - // initializes the point source id counter - this->clear(); - bulkSourceIds_.clear(); - bulkSourceWeights_.clear(); - extendedSourceStencil_.clear(); - - // precompute the vertex indices for efficiency for the box method - this->precomputeVertexIndices(bulkIdx); - this->precomputeVertexIndices(lowDimIdx); + // prepare the internal data structures + prepareDataStructures_(); const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); - bulkSourceIds_.resize(this->gridView(bulkIdx).size(0)); - bulkSourceWeights_.resize(this->gridView(bulkIdx).size(0)); - - // intersect the bounding box trees - this->glueGrids(); - - // reserve memory for data - this->pointSourceData().reserve(this->glue().size()); - this->averageDistanceToBulkCell().reserve(this->glue().size()); - fluxScalingFactor_.reserve(this->glue().size()); - - // reserve memory for stencils - const auto numBulkElements = this->gridView(bulkIdx).size(0); - for (GridIndex bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) - { - this->couplingStencils(bulkIdx)[bulkElementIdx].reserve(50); - extendedSourceStencil_.stencil()[bulkElementIdx].reserve(50); - bulkSourceIds_[bulkElementIdx][0].reserve(20); - bulkSourceWeights_[bulkElementIdx][0].reserve(20); - } - // generate a bunch of random vectors and values for // Monte-carlo integration on the cylinder defined by line and radius static const int samples = getParam("MixedDimension.NumMCSamples"); @@ -1301,39 +1272,8 @@ public: } } - // make extra stencils unique - for (auto&& stencil : extendedSourceStencil_.stencil()) - { - std::sort(stencil.second.begin(), stencil.second.end()); - stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); - - // remove the vertices element (box) - if (isBox()) - { - const auto& indices = this->vertexIndices(bulkIdx, stencil.first); - stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), - [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), - stencil.second.end()); - } - // remove the own element (cell-centered) - else - { - stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), - [&](auto i){ return i == stencil.first; }), - stencil.second.end()); - } - } - - // make stencils unique - using namespace Dune::Hybrid; - forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) - { - for (auto&& stencil : this->couplingStencils(domainIdx)) - { - std::sort(stencil.second.begin(), stencil.second.end()); - stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); - } - }); + // make the stencils unique + makeUniqueStencil_(); if (!this->pointSources(bulkIdx).empty()) DUNE_THROW(Dune::InvalidStateException, "Kernel method shouldn't have point sources in the bulk domain but only volume sources!"); @@ -1471,6 +1411,79 @@ private: return integral/length; } + void prepareDataStructures_() + { + // clear all internal members like pointsource vectors and stencils + // initializes the point source id counter + this->clear(); + bulkSourceIds_.clear(); + bulkSourceWeights_.clear(); + extendedSourceStencil_.stencil().clear(); + + // precompute the vertex indices for efficiency for the box method + this->preComputeVertexIndices(bulkIdx); + this->preComputeVertexIndices(lowDimIdx); + + bulkSourceIds_.resize(this->gridView(bulkIdx).size(0)); + bulkSourceWeights_.resize(this->gridView(bulkIdx).size(0)); + + // intersect the bounding box trees + this->glueGrids(); + + // reserve memory for data + this->pointSourceData().reserve(this->glue().size()); + this->averageDistanceToBulkCell().reserve(this->glue().size()); + fluxScalingFactor_.reserve(this->glue().size()); + + // reserve memory for stencils + const auto numBulkElements = this->gridView(bulkIdx).size(0); + for (GridIndex bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) + { + this->couplingStencils(bulkIdx)[bulkElementIdx].reserve(50); + extendedSourceStencil_.stencil()[bulkElementIdx].reserve(50); + bulkSourceIds_[bulkElementIdx][0].reserve(20); + bulkSourceWeights_[bulkElementIdx][0].reserve(20); + } + } + + //! Make the stencils unique + void makeUniqueStencil_() + { + // make extra stencils unique + for (auto&& stencil : extendedSourceStencil_.stencil()) + { + std::sort(stencil.second.begin(), stencil.second.end()); + stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + + // remove the vertices element (box) + if (isBox()) + { + const auto& indices = this->vertexIndices(bulkIdx, stencil.first); + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), + stencil.second.end()); + } + // remove the own element (cell-centered) + else + { + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return i == stencil.first; }), + stencil.second.end()); + } + } + + // make stencils unique + using namespace Dune::Hybrid; + forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) + { + for (auto&& stencil : this->couplingStencils(domainIdx)) + { + std::sort(stencil.second.begin(), stencil.second.end()); + stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + } + }); + } + //! add additional stencil entries for the bulk element void addBulkSourceStencils_(GridIndex bulkElementIdx, GridIndex coupledLowDimElementIdx, GridIndex coupledBulkElementIdx) { -- GitLab From 778569671649551546a86de1c86e40c008d46847 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Mon, 15 Jul 2019 18:51:44 +0200 Subject: [PATCH 16/17] Add embedded coupling modes --- dumux/multidomain/embedded/couplingmanager1d3d.hh | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index a477bdba3a..c6d09501f6 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -51,7 +51,7 @@ namespace Dumux { * \brief The coupling mode */ enum class EmbeddedCouplingMode -{ line, average, cylindersources, kernel, projection }; +{ line, average, cylindersources, kernel, kernelAnisotropic, wellIndex, projection }; /*! * \ingroup EmbeddedCoupling @@ -65,6 +65,8 @@ std::string toString(EmbeddedCouplingMode mode) case EmbeddedCouplingMode::average: return "average"; case EmbeddedCouplingMode::cylindersources: return "cylinder"; case EmbeddedCouplingMode::kernel: return "kernel"; + case EmbeddedCouplingMode::kernelAnisotropic: return "kernelAnisotropic"; + case EmbeddedCouplingMode::wellIndex: return "wellIndex"; case EmbeddedCouplingMode::projection: return "projection"; default: return "unknown"; } -- GitLab From 7346d9398165dbabca8ea687a836019b76dfd969 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Sat, 14 Dec 2019 15:32:20 +0100 Subject: [PATCH 17/17] [1d3d] Optimize cylinder source coupling mode --- .../embedded/couplingmanager1d3d.hh | 75 ++++++++++--------- 1 file changed, 41 insertions(+), 34 deletions(-) diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index c6d09501f6..316047d549 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -28,6 +28,7 @@ #include #include +#include #include #include @@ -390,15 +391,14 @@ public: // precompute interpolation data for box scheme for each cut bulk element if constexpr (isBox()) { - if (!static_cast(circleCornerIndices.count(bulkElementIdx))) - { - const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); - circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx); + const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); + circleCornerIndices.push_back(&(this->vertexIndices(bulkIdx, bulkElementIdx))); - // evaluate shape functions at the integration point - const auto bulkGeometry = bulkElement.geometry(); - this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); - } + // evaluate shape functions at the integration point + const auto bulkGeometry = bulkElement.geometry(); + ShapeValues shapeValues; + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], shapeValues); + circleShapeValues.emplace_back(std::move(shapeValues)); } } } @@ -721,6 +721,9 @@ public: */ void computePointSourceData(std::size_t order = 1, bool verbose = false) { + // if we use the circle average as the 3D values or a point evaluation + static const bool useCircleAverage = getParam("MixedDimension.UseCircleAverage", true); + // Initialize the bulk bounding box tree const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); @@ -792,11 +795,11 @@ public: std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); std::vector> circleStencil; circleStencil.reserve(circlePoints.size()); // for box - std::unordered_map, std::vector> > circleCornerIndices; + std::vector>*> circleCornerIndices; using ShapeValues = std::vector >; - std::unordered_map, ShapeValues> circleShapeValues; + std::vector circleShapeValues; - // go over all circle points and precompute some quantities for the circle average + // go over all points of the average operator for (int k = 0; k < circlePoints.size(); ++k) { circleBulkElementIndices[k] = intersectingEntities(circlePoints[k], bulkTree); @@ -809,18 +812,17 @@ public: circleStencil.push_back(bulkElementIdx); circleIpWeight.push_back(localCircleAvgWeight); - // precompute interpolation data for the box scheme + // precompute interpolation data for box scheme for each cut bulk element if (isBox()) { - if (!static_cast(circleCornerIndices.count(bulkElementIdx))) - { - const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); - circleCornerIndices[bulkElementIdx] = this->vertexIndices(bulkIdx, bulkElementIdx); + const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); + circleCornerIndices.push_back(&(this->vertexIndices(bulkIdx, bulkElementIdx))); - // evaluate shape functions at the integration point - const auto bulkGeometry = bulkElement.geometry(); - this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], circleShapeValues[bulkElementIdx]); - } + // evaluate shape functions at the integration point + const auto bulkGeometry = bulkElement.geometry(); + ShapeValues shapeValues; + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], shapeValues); + circleShapeValues.emplace_back(std::move(shapeValues)); } } } @@ -832,7 +834,7 @@ public: for (const auto& vertices : circleCornerIndices) { this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), - vertices.second.begin(), vertices.second.end()); + vertices->begin(), vertices->end()); } } @@ -880,7 +882,8 @@ public: // add data needed to compute integral over the circle if constexpr (isBox()) { - psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil); + if (useCircleAverage) + psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil); using ShapeValues = std::vector >; const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); @@ -890,7 +893,8 @@ public: } else { - psData.addCircleInterpolation(circleIpWeight, circleStencil); + if (useCircleAverage) + psData.addCircleInterpolation(circleIpWeight, circleStencil); psData.addBulkInterpolation(bulkElementIdx); } @@ -911,21 +915,24 @@ public: this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); } - // export bulk circle stencil - if (isBox()) + // export bulk circle stencil (only needed for circle average) + if (useCircleAverage) { - // we insert all vertices and make it unique later - for (const auto& vertices : circleCornerIndices) + if (isBox()) { - extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), - vertices.second.begin(), vertices.second.end()); + // we insert all vertices and make it unique later + for (const auto& vertices : circleCornerIndices) + { + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + vertices->begin(), vertices->end()); + } + } + else + { + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + circleStencil.begin(), circleStencil.end()); } - } - else - { - extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), - circleStencil.begin(), circleStencil.end()); } } } -- GitLab