Additional file needed for the CO2 test in order to get intersections

from vertices.


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9132 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/boxmodels/common/intersectiontovertexbc.hh

+/*****************************************************************************
+ *   Copyright (C) 2010 by Bernd Flemisch                                    *
+ *   Institute for Modelling Hydraulic and Environmental Systems             *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for all problems which use the box scheme
+ */
+#ifndef DUMUX_INTERSECTIONTOVERTEXBC_HH
+#define DUMUX_INTERSECTIONTOVERTEXBC_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+
+namespace Dumux
+{
+
+template<typename TypeTag> 
+class IntersectionToVertexBC
+{
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+    enum {
+        numEq = GET_PROP_VALUE(TypeTag, NumEq),
+        dim = GridView::dimension,
+    };
+
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+    typedef typename GridView::IntersectionIterator IntersectionIterator;
+    typedef typename GridView::template Codim<dim>::Entity Vertex;
+
+    typedef typename Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+    typedef typename Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
+
+public:
+    IntersectionToVertexBC(const Problem& problem)
+    : problem_(problem)
+    {
+	vertexBC.resize(problem.vertexMapper().size());
+        for (int vertIdx = 0; vertIdx < vertexBC.size(); vertIdx++)
+            vertexBC[vertIdx].setAllNeumann();
+
+        ElementIterator eIt = problem.gridView().template begin<0>();
+        const ElementIterator elemEndIt = problem.gridView().template end<0>();
+        for (; eIt != elemEndIt; ++eIt) {
+            Dune::GeometryType geoType = eIt->geometry().type();
+            const ReferenceElement &refElem = ReferenceElements::general(geoType);
+
+            IntersectionIterator isIt = problem.gridView().ibegin(*eIt);
+            IntersectionIterator isEndIt = problem.gridView().iend(*eIt);
+            for (; isIt != isEndIt; ++isIt) {
+                if (!isIt->boundary())
+                    continue;
+
+                BoundaryTypes bcTypes;
+		problem.boundaryTypes(bcTypes, *isIt);
+
+		if (!bcTypes.hasDirichlet())
+                    continue;
+
+                int faceIdx = isIt->indexInInside();
+                int numFaceVerts = refElem.size(faceIdx, 1, dim);
+                for (int faceVertIdx = 0; faceVertIdx < numFaceVerts; ++faceVertIdx)
+                {
+                    int elemVertIdx = refElem.subEntity(faceIdx, 1, faceVertIdx, dim);
+                    int globalVertIdx = problem.vertexMapper().map(*eIt, elemVertIdx, dim);
+
+                    for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+                      if (bcTypes.isDirichlet(eqIdx))
+                          vertexBC[globalVertIdx].setDirichlet(eqIdx);
+                }
+            }
+        }
+    }
+
+    void boundaryTypes(BoundaryTypes& values, const Vertex& vertex) const
+    {
+        values.setAllNeumann();
+
+	int vertIdx = problem_.vertexMapper().map(vertex);
+	const BoundaryTypes& bcTypes = vertexBC[vertIdx];
+
+        for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+            if (bcTypes.isDirichlet(eqIdx))
+                values.setDirichlet(eqIdx);
+    }
+
+private: 
+    const Problem& problem_;
+    std::vector<BoundaryTypes> vertexBC;
+};
+}
+
+#endif
+