diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 83bbba01e6b484955562fea2179f526753f21739..b2ce43e3b6414bafaf87d64d132478401d8ed2b6 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -348,6 +348,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
{
entries = 0.;
ElementPtr elementPtrI = intersection.inside();
+ int globalIdxI = problem().variables().index(*elementPtrI);
// get global coordinate of cell center
const GlobalPosition& globalPos = elementPtrI->geometry().center();
@@ -509,8 +510,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
{
if (cellDataI.wasRefined() && cellDataJ.wasRefined())
{
- problem().variables().cellData(problem().variables().index(*elementPtrI)).setUpwindCell(intersection.indexInInside(),
- contiWEqIdx, false);
+ problem().variables().cellData(globalIdxI).setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
}
@@ -543,8 +543,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
{
if (cellDataI.wasRefined() && cellDataJ.wasRefined())
{
- problem().variables().cellData(problem().variables().index(*elementPtrI)).setUpwindCell(intersection.indexInInside(),
- contiNEqIdx, false);
+ problem().variables().cellData(globalIdxI).setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
}
Scalar averagedMassFraction[2];
diff --git a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
index e2b55e8a3f72c917379fbdc2cf9388d55d14da52..d16ec39cbdb1269f002b4ab0e52c0d5fb56d9e43 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
@@ -359,8 +359,8 @@ void FVPressure2P2CMultiPhysics<TypeTag>::assemble(bool first)
* \param cellDataI Data of cell I
*/
template<class TypeTag>
-void FVPressure2P2CMultiPhysics<TypeTag>::get1pSource(Dune::FieldVector<Scalar,
- 2>& sourceEntry, const Element& elementI, const CellData& cellDataI)
+void FVPressure2P2CMultiPhysics<TypeTag>::get1pSource(Dune::FieldVector<Scalar, 2>& sourceEntry,
+ const Element& elementI, const CellData& cellDataI)
{
sourceEntry=0.;
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index f884e49ba29b11857e504311d5258f93fea2b1cd..7271d78a85373fcb9ecd04f530f4290d9ce8ac93 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -284,7 +284,8 @@ public:
localTimeStepping_ = subCFLFactor_/cFLFactor < 1.0 - dtThreshold_;
if (localTimeStepping_)
- std::cout<<"max CFL-Number of "<<cFLFactor<<", max Sub-CFL-Number of "<<subCFLFactor_<<": Enable local time-stepping!\n";
+ std::cout<<"max CFL-Number of "<<cFLFactor<<", max Sub-CFL-Number of "
+ <<subCFLFactor_<<": Enable local time-stepping!" << std::endl;
}
virtual ~FVTransport2P2C()
@@ -768,8 +769,8 @@ void FVTransport2P2C<TypeTag>::getFlux(ComponentVector& fluxEntries,
// verbose (only for one side)
if(globalIdxI > globalIdxJ)
Dune::dinfo << "harmonicMean flux of phase" << phaseIdx <<" used from cell" << globalIdxI<< " into " << globalIdxJ
- << " ; TE upwind I = "<< cellDataI.isUpwindCell(intersection.indexInInside(),
- contiEqIdx) << " but pot = "<< potential[phaseIdx] << " \n";
+ << " ; TE upwind I = "<< cellDataI.isUpwindCell(intersection.indexInInside(), contiEqIdx)
+ << " but pot = "<< potential[phaseIdx] << std::endl;
#endif
}
@@ -1330,7 +1331,7 @@ void FVTransport2P2C<TypeTag>::innerUpdate(TransportSolutionType& updateVec)
subDt += dtCorrection;
if (verbosity_ > 0)
- std::cout<<" Sub-time-step size: "<<subDt<<"\n";
+ std::cout<<" Sub-time-step size: "<<subDt<< std::endl;
bool stopTimeStep = false;
int size = problem_.gridView().size(0);
diff --git a/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh b/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh
index ad8818407b7e747edee11e92d454a6d83549049a..1f0a77c9528af73ed4172633dcb007e4fdd50748 100644
--- a/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh
+++ b/dumux/decoupled/2p2c/variableclass2p2cadaptive.hh
@@ -37,9 +37,9 @@ namespace Dumux
//! Class holding additionally mpfa data of adaptive compositional models.
/*!
* This class provides the possibility to store and load the transmissibilities (and associated infos)
- * of the mpfa method per irregular face. While this class provides the storage container for both 2D
- * and 3D implementation, only the 2D storage methods are included here.
- * Note that according to the number of half-edges (sub-faces) regarded, one ore two transmissibility
+ * of the mpfa method per irregular face. This class provides the storage container and access methods
+ * for both 2D and 3D implementation.
+ * Note that according to the number of half-edges (sub-faces) regarded, one ore more transmissibility
* can be stored and loaded.
*
* @tparam TypeTag The Type Tag
@@ -67,7 +67,10 @@ private:
};
enum //!< for first and second half edge (2D) or subvolume face (3D)
{
- first = 0, second = 1
+ first = 0,
+ second = 1,
+ diagonal1 = 2,
+ diagonal2 = 3
};
// convenience shortcuts for Vectors/Matrices
typedef Dune::FieldVector<Scalar, GridView::dimensionworld> GlobalPosition;
@@ -75,40 +78,40 @@ private:
protected:
/** in the 2D case, we need to store 1 additional cell. In 3D, we store
- * dim-1=2 cells for one interaction volume, and 4 if two interaction volumes
+ * dim-1=2 cells for one interaction volume, and 8 cells if four interaction volumes
* are regarded.
*/
- const static int storageRequirement = (dim-1)*(dim-1);
+ const static int storageRequirement = pow((dim-1), dim);
//! Storage object for data related to the MPFA method
/**
- * This Struct stores the transmissibility coefficients
- * for the two half-eges of an irregular interface (one
- * near a hanging node) in an h-adaptive simulation.
+ * This Object stores the transmissibility coefficients
+ * for all interaction regions/volumes of an irregular interface
+ * in an h-adaptive simulation. It is valid for 2D/3D
*/
struct mpfaData
{
- TransmissivityMatrix T1_[2];
+ TransmissivityMatrix T1_[(dim-1)*2];
int globalIdx3_[storageRequirement];
GlobalPosition globalPos3_[storageRequirement];
std::vector<IntersectionIterator> secondHalfEdgeIntersection_;
- bool hasSecondHalfEdge;
+ int interactionRegionsStored;
//! Constructor for the local storage object of mpfa data
mpfaData()
{
- hasSecondHalfEdge = false;
+ interactionRegionsStored = 0;
}
- //! stores an intersection
- /** This also provides the information that both half-edges are
- * regarded and information was stored.
+ //! stores an intersection for the 2D implementation
+ /** This method also provides the information that both half-edges (2D) are
+ * regarded and information was stored: Two interaction regions are applied
* \param is23 Intersection pointing to 3rd cell of additional interaction region
*/
void setIntersection(IntersectionIterator& is23)
{
secondHalfEdgeIntersection_.push_back(is23);
- hasSecondHalfEdge = true;
+ interactionRegionsStored = 2;
};
- //! Acess method to the stored intersection
+ //! Acess method to the stored intersection (only 2D)
const IntersectionIterator& getIntersection()
{
return secondHalfEdgeIntersection_[0];
@@ -153,8 +156,8 @@ public:
problem.pressureModel().adaptPressure();
}
- //! Stores Mpfa Data on an intersection
- /** The method stores information to the interaction region (Transmissitivity
+ //! Stores Mpfa Data of one interaction region on an intersection
+ /** The method stores information of ONE interaction region (Transmissitivity
* as well as details of the 3rd cell in the region) into a storage container.
* The key to each element is the index of the intersection, seen from the smaller
* cell (only this is unique).
@@ -192,8 +195,9 @@ public:
irregularInterfaceMap_[intersectionID].globalPos3_[0] = globalPos3;
irregularInterfaceMap_[intersectionID].globalIdx3_[0] = globalIdx3;
- return;
}
+
+
//! Stores Mpfa Data on an intersection for both half-edges
/** The method stores information of both interaction regions (Transmissitivity
* as well as details of the 3rd cells of both regions) into a storage container.
@@ -253,9 +257,110 @@ public:
return;
}
- //! Provides acess to stored Mpfa Data on an intersection for both half-edges
+ //! Stores 3D Mpfa Data on an intersection
+ /** The method stores information of the interaction region (Transmissitivity
+ * as well as details of the 3rd & 4th cell in the region) into a storage container.
+ * The key to each element is the index of the intersection, seen from the smaller
+ * cell (only this index is unique).
+ * If we arrive from the "wrong" (i.e. non-unique) direction, we invert fluxes.
+ *
+ * \param irregularIs The current irregular intersection
+ * \param T1 Transmissitivity matrix for flux calculations
+ * \param globalPos3 The position of the 3rd cell of the interaction region
+ * \param globalIdx3 The index of the 3rd cell of the interaction region
+ * \param globalPos4 The position of the 4th cell of the interaction region
+ * \param globalIdx4 The index of the 4th cell of the interaction region
+ * \param subFaceIdx The index of the subface (up to 4 subfaces possible in 3D)
+ */
+ void storeMpfaData3D(const typename GridView::Intersection & irregularIs,
+ const TransmissivityMatrix& T1,
+ const GlobalPosition& globalPos3,
+ const int& globalIdx3,
+ const GlobalPosition& globalPos4,
+ const int& globalIdx4,
+ int subFaceIdx = 0)
+ {
+ // if the second interaction Region (subfaceIdx=1) needs to be stored,
+ // we need an offset for the second globalPos and globalIdxes
+ const int offset = subFaceIdx * 2;
+
+ IdType intersectionID
+ = grid_.localIdSet().subId(*irregularIs.inside(),
+ irregularIs.indexInInside(), 1);
+
+ // mapping is only unique from smaller cell (if *inside and not *outside)
+ if (irregularIs.inside().level() < irregularIs.outside().level())
+ {
+ // IS is regarded from larger cell: get the unique ID as seen from smaller
+ intersectionID
+ = grid_.localIdSet().subId(*irregularIs.outside(),
+ irregularIs.indexInOutside(), 1);
+
+ // store as if it was seen from smaller: change i & j
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][0] = -T1[1];
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][1] = -T1[0];
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][2] = -T1[2];
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][3] = -T1[3];
+ }
+ else
+ {
+ // proceed with numbering according to case2
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][0] = T1[0];
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][1] = T1[1];
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][2] = T1[2];
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][3] = T1[3];
+ }
+ irregularInterfaceMap_[intersectionID].globalPos3_[offset+0] = globalPos3;
+ irregularInterfaceMap_[intersectionID].globalIdx3_[offset+0] = globalIdx3;
+ irregularInterfaceMap_[intersectionID].globalPos3_[offset+1] = globalPos4;
+ irregularInterfaceMap_[intersectionID].globalIdx3_[offset+1] = globalIdx4;
+
+ irregularInterfaceMap_[intersectionID].interactionRegionsStored
+ = std::max(irregularInterfaceMap_[intersectionID].interactionRegionsStored, subFaceIdx+1);
+
+ return;
+ }
+
+ //! Weigths the transmissivity coefficient by the flux area (3D)
+ /** Each cell face could contain up to 4 interaction regions in 3D. If fewer interaction
+ * regions are considered (e.g. at a boundary), the flux goes through a higher area than
+ * was regarded in the calculation of the interaction region. Therefore, the transmissibility
+ * coefficients have to be weighted (scaled).
+ * If no specific face Index is given, all present Transmissibility matrices are weighted
+ * homogeneously by the given weight.
+ * \param irregularIs The current irregular intersection
+ * \param weight The weighting factor
+ * \param subFaceIdx The index of the subface (up to 4 subfaces possible in 3D)
+ */
+ void performTransmissitivityWeighting(const typename GridView::Intersection & irregularIs,
+ Scalar weight, int subFaceIdx = -1)
+ {
+ IdType intersectionID
+ = grid_.localIdSet().subId(*irregularIs.inside(),
+ irregularIs.indexInInside(), 1);
+
+ // mapping is only unique from smaller cell (if *inside and not *outside)
+ if (irregularIs.inside().level() < irregularIs.outside().level())
+ {
+ // IS is regarded from larger cell: get the unique ID as seen from smaller
+ intersectionID
+ = grid_.localIdSet().subId(*irregularIs.outside(),
+ irregularIs.indexInOutside(), 1);
+ }
+
+ // for subFaceIdx == -1, we weight all subfaces equally
+ if(subFaceIdx == -1)
+ {
+ for(int i = 0; i < irregularInterfaceMap_[intersectionID].interactionRegionsStored; i++)
+ irregularInterfaceMap_[intersectionID].T1_[i] *= weight;
+ }
+ else //weight specifically
+ irregularInterfaceMap_[intersectionID].T1_[subFaceIdx] *= weight;
+ }
+
+ //! Provides access to stored Mpfa Data on an intersection for both half-edges
/** The method gets information of both interaction regions (Transmissitivity
- * as well as details of the 3rd cells of both regions) into a storage container.
+ * as well as details of the 3rd cells of both regions) from a storage container.
* The key to each element is the index of the intersection, seen from the smaller
* cell (only this is unique).
* If we arrive from the "wrong" (i.e. non-unique) direction, we invert fluxes.
@@ -298,7 +403,7 @@ public:
globalPos3 = irregularInterfaceMap_[intersectionID].globalPos3_[0];
globalIdx3 = irregularInterfaceMap_[intersectionID].globalIdx3_[0];
//second half edge
- if(irregularInterfaceMap_[intersectionID].hasSecondHalfEdge)
+ if(irregularInterfaceMap_[intersectionID].interactionRegionsStored == 2)
{
secondHalfEdgeIntersectionIt = irregularInterfaceMap_[intersectionID].getIntersection();
T1_secondHalfEdge[0] = -irregularInterfaceMap_[intersectionID].T1_[second][2];
@@ -320,7 +425,7 @@ public:
globalPos3 = irregularInterfaceMap_[intersectionID].globalPos3_[0];
globalIdx3 = irregularInterfaceMap_[intersectionID].globalIdx3_[0];
//second half edge
- if(irregularInterfaceMap_[intersectionID].hasSecondHalfEdge)
+ if(irregularInterfaceMap_[intersectionID].interactionRegionsStored == 2)
{
secondHalfEdgeIntersectionIt = irregularInterfaceMap_[intersectionID].getIntersection();
T1_secondHalfEdge = irregularInterfaceMap_[intersectionID].T1_[second];
@@ -331,6 +436,81 @@ public:
}
return 1;
}
+
+
+ //! Provides access to stored 3D Mpfa Data on an intersection for up to 4 subfaces
+ /** The method gets information of up to 4 interaction regions (Transmissitivity
+ * as well as details of the 3rd & 4th cells of the regions) from a storage container.
+ * The key to each element is the index of the intersection, seen from the smaller
+ * cell (only this is unique).
+ * If we arrive from the "wrong" (i.e. non-unique) direction, we invert fluxes.
+ *
+ * \param irregularIs The current irregular intersection
+ * \param secondHalfEdgeIntersectionIt Iterator to the intersection connecting the second interaction region
+ * \param T1 Transmissitivity matrix for flux calculations: unique interaction region
+ * \param T1_secondHalfEdge Second transmissitivity matrix for flux calculations for non-unique interaction region
+ * \param globalPos3 The position of the 3rd cell of the first interaction region
+ * \param globalIdx3 The index of the 3rd cell of the first interaction region
+ * \param globalPos4 The position of the 4th cell of the interaction region
+ * \param globalIdx4 The index of the 4th cell of the interaction region
+ * \param subFaceIdx The index of the subface (up to 4 subfaces possible in 3D)
+ */
+ int getMpfaData3D(const Intersection& irregularIs,
+ TransmissivityMatrix& T1,
+ GlobalPosition& globalPos3,
+ int& globalIdx3,
+ GlobalPosition& globalPos4,
+ int& globalIdx4,
+ int subFaceIdx = 0)
+ {
+ // if the second interaction Region (subfaceIdx=1) needs to be stored,
+ // we need an offset for the second globalPos and globalIdxes
+ const int offset = subFaceIdx * 2;
+
+ IdType intersectionID
+ = grid_.localIdSet().subId(*irregularIs.inside(),
+ irregularIs.indexInInside(), 1);
+ // mapping is only unique from smaller cell (if *inside and not *outside)
+ if (irregularIs.inside().level() < irregularIs.outside().level())
+ {
+ // IS is regarded from larger cell: get the unique number as seen from smaller
+ intersectionID
+ = grid_.localIdSet().subId(*irregularIs.outside(),
+ irregularIs.indexInOutside(), 1);
+
+ // check if T1ransmissibility matrix was stored for that IF
+ if (irregularInterfaceMap_.find(intersectionID) == irregularInterfaceMap_.end())
+ return 0; // no stored data!
+
+ // If data is stored, it is so as if the IF is regarded from smaller cell.
+ // since we are looking from larger cell, cells i & j have to be changed
+ // Additionally, flux points in opposite direction: - sign
+ T1[0] = -irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][1];
+ T1[1] = -irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][0];
+ T1[2] = -irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][2];
+ T1[3] = -irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][3];
+
+ }
+ else
+ {
+ // check if T1ransmissibility matrix was stored for that IF
+ if (irregularInterfaceMap_.find(intersectionID) == irregularInterfaceMap_.end())
+ return 0; // no stored data!
+
+ T1[0] = irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][0];
+ T1[1] = irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][1];
+ T1[2] = irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][2];
+ T1[3] = irregularInterfaceMap_[intersectionID].T1_[subFaceIdx][3];
+ }
+
+ // return what does not depend on direction: additional cells
+ globalPos3 = irregularInterfaceMap_[intersectionID].globalPos3_[offset+0];
+ globalIdx3 = irregularInterfaceMap_[intersectionID].globalIdx3_[offset+0];
+ globalPos4 = irregularInterfaceMap_[intersectionID].globalPos3_[offset+1];
+ globalIdx4 = irregularInterfaceMap_[intersectionID].globalIdx3_[offset+1];
+
+ return irregularInterfaceMap_[intersectionID].interactionRegionsStored;
+ }
//@}
};