test_gridmanager.cc 14.9 KB
Newer Older
1
2
3
4
5
6
7
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   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    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
 *   (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
21
 * \ingroup FacetTests
22
 * \brief Tests the grid manager class for models using facet coupling.
23
 */
24

25
26
27
28
29
30
31
#include <config.h>
#include <iostream>

#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>

#include <dune/grid/uggrid.hh>
32
#include <dune/grid/common/mcmgmapper.hh>
33
34
35
36
37
38
#include <dune/alugrid/grid.hh>
#include <dune/foamgrid/foamgrid.hh>

#include <dune/grid/io/file/vtk/vtkwriter.hh>

#include <dumux/common/parameters.hh>
39
#include <dumux/io/vtk/function.hh>
40
#include <dumux/multidomain/facet/gridmanager.hh>
41
42
43
44
45

#ifndef BULKGRIDTYPE // default to ug grid if not provided by CMake
#define BULKGRIDTYPE Dune::UGGrid<3>
#endif

46
int main (int argc, char *argv[])
47
48
49
50
51
52
53
{
    // maybe initialize mpi
    Dune::MPIHelper::instance(argc, argv);

    // parse command line argument parameters
    Dumux::Parameters::init(argc, argv);

54
    // make grid from input file
55
56
57
    using BulkGrid = BULKGRIDTYPE;
    using FacetGrid = Dune::FoamGrid<2, 3>;
    using EdgeGrid = Dune::FoamGrid<1, 3>;
58
59
60
    using GridManager = Dumux::FacetCouplingGridManager<BulkGrid, FacetGrid, EdgeGrid>;
    GridManager gridManager;
    gridManager.init();
61
62

    // check grid sizes
63
64
65
    const auto& bulkGridView = gridManager.grid<0>().leafGridView();
    const auto& facetGridView = gridManager.grid<1>().leafGridView();
    const auto& edgeGridView = gridManager.grid<2>().leafGridView();
66

67
68
    if (bulkGridView.size(0) != 491)
        DUNE_THROW(Dune::InvalidStateException, "Bulk grid has " << bulkGridView.size(0) << " instead of 491 elements");
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
    if (bulkGridView.size(3) != 153)
        DUNE_THROW(Dune::InvalidStateException, "Bulk grid has " << bulkGridView.size(3) << " instead of 153 vertices");

    if (facetGridView.size(0) != 32)
        DUNE_THROW(Dune::InvalidStateException, "Facet grid has " << facetGridView.size(0) << " instead of 32 elements");
    if (facetGridView.size(2) != 23)
        DUNE_THROW(Dune::InvalidStateException, "Facet grid has " << facetGridView.size(2) << " instead of 23 vertices");

    if (edgeGridView.size(0) != 2)
        DUNE_THROW(Dune::InvalidStateException, "Edge grid has " << edgeGridView.size(0) << " instead of 2 elements");
    if (edgeGridView.size(1) != 3)
        DUNE_THROW(Dune::InvalidStateException, "Edge grid has " << edgeGridView.size(1) << " instead of 3 vertices");

    // all entities should have the domain index 10
    for (const auto& e : elements(bulkGridView))
    {
85
        const auto domainMarker = gridManager.getGridData()->getElementDomainMarker<0>(e);
86
87
88
89
90
91
        if (domainMarker != 10)
            DUNE_THROW(Dune::InvalidStateException, "Bulk element has domain marker " << domainMarker << " instead of 10");
    }

    for (const auto& e : elements(facetGridView))
    {
92
        const auto domainMarker = gridManager.getGridData()->getElementDomainMarker<1>(e);
93
94
95
96
97
98
        if (domainMarker != 10)
            DUNE_THROW(Dune::InvalidStateException, "Facet element has domain marker " << domainMarker << " instead of 10");
    }

    for (const auto& e : elements(edgeGridView))
    {
99
        const auto domainMarker = gridManager.getGridData()->getElementDomainMarker<2>(e);
100
101
102
103
        if (domainMarker != 10)
            DUNE_THROW(Dune::InvalidStateException, "Edge element has domain marker " << domainMarker << " instead of 10");
    }

Hanchuan Wu's avatar
Hanchuan Wu committed
104
105
    // Check boundary segments. There should be mesh file-defined
    // segments on all sides except the positive x-direction
106
    std::size_t bSegmentCount = 0;
107
108
    std::size_t posXCount, negXCount, posYCount, negYCount, posZCount, negZCount;
    posXCount = negXCount = posYCount = negYCount = posZCount = negZCount = 0;
109
110
111
112
113
114
115
116
    for (const auto& e : elements(bulkGridView))
    {
        for (const auto& is : intersections(bulkGridView, e))
        {
            if (is.boundary())
            {
                const auto n = is.centerUnitOuterNormal();

117
                if (gridManager.getGridData()->wasInserted<0>(is))
118
119
                {
                    bSegmentCount++;
120
                    const auto insIdx = gridManager.getEmbeddings()->insertionIndex<0>(is);
121
122
123
                    const auto marker = gridManager.getGridData()->getBoundaryDomainMarker<0>(insIdx);

                    if (Dune::FloatCmp::eq(n[0], 1.0, 1e-6)) // pos x-dir
124
                    {
125
126
127
128
129
130
131
                        if (marker != 1)
                            DUNE_THROW(Dune::InvalidStateException, "positive x-face should have marker 7, but has " << marker);
                        posXCount++;
                    }
                    else if (Dune::FloatCmp::eq(n[0], -1.0, 1e-6)) // neg x-dir
                    {
                        if (marker != 2)
132
133
134
135
136
                            DUNE_THROW(Dune::InvalidStateException, "negative x-face should have marker 4, but has " << marker);
                        negXCount++;
                    }
                    else if (Dune::FloatCmp::eq(n[1], 1.0, 1e-6)) // pos y-dir
                    {
137
                        if (marker != 3)
138
139
140
141
142
                            DUNE_THROW(Dune::InvalidStateException, "positive y-face should have marker 5, but has " << marker);
                        posYCount++;
                    }
                    else if (Dune::FloatCmp::eq(n[1], -1.0, 1e-6)) // neg y-dir
                    {
143
                        if (marker != 4)
144
145
146
147
148
                            DUNE_THROW(Dune::InvalidStateException, "negative y-face should have marker 6, but has " << marker);
                        negYCount++;
                    }
                    else if (Dune::FloatCmp::eq(n[2], 1.0, 1e-6)) // pos z-dir
                    {
149
                        if (marker != 5)
150
151
152
153
154
                            DUNE_THROW(Dune::InvalidStateException, "positive z-face should have marker 1, but has " << marker);
                        posZCount++;
                    }
                    else if (Dune::FloatCmp::eq(n[2], -1.0, 1e-6)) // neg z-dir
                    {
155
                        if (marker != 6)
156
157
158
                            DUNE_THROW(Dune::InvalidStateException, "negative z-face should have marker 2, but has " << marker);
                        negZCount++;
                    }
159
160
                    else
                        DUNE_THROW(Dune::InvalidStateException, "Could not deduce boundary segment orientation");
161
162
                }
                else
163
                    DUNE_THROW(Dune::InvalidStateException, "Boundary intersection was not inserted. Can't obtain boundary segment index");
164
165
166
167
168
            }
        }
    }

    //! make sure we found the right number of boundary segments
169
170
    if (bSegmentCount != 240) DUNE_THROW(Dune::InvalidStateException, "Found " << bSegmentCount << " instead of 200 boundary segments");
    if (posXCount != 40) DUNE_THROW(Dune::InvalidStateException, "Found " << posXCount << " instead of 40 boundary segments in pos x-direction");
171
172
173
174
175
176
177
    if (negXCount != 40) DUNE_THROW(Dune::InvalidStateException, "Found " << negXCount << " instead of 40 boundary segments in neg x-direction");
    if (posYCount != 40) DUNE_THROW(Dune::InvalidStateException, "Found " << posYCount << " instead of 40 boundary segments in pos y-direction");
    if (negYCount != 40) DUNE_THROW(Dune::InvalidStateException, "Found " << negYCount << " instead of 40 boundary segments in neg y-direction");
    if (posZCount != 40) DUNE_THROW(Dune::InvalidStateException, "Found " << posZCount << " instead of 40 boundary segments in pos z-direction");
    if (negZCount != 40) DUNE_THROW(Dune::InvalidStateException, "Found " << negZCount << " instead of 40 boundary segments in neg z-direction");

    //! check if we found the right number of embeddings
178
179
    const auto& edgeEmbeddings = gridManager.getEmbeddings()->embeddedEntityMap(2);
    const auto& edgeEmbedments = gridManager.getEmbeddings()->adjoinedEntityMap(2);
180
181
182
183
184
185
    if (edgeEmbeddings.size() != 0) DUNE_THROW(Dune::InvalidStateException, "The grid with lowest dimension can't have embedded entities");
    if (edgeEmbedments.size() != 2) DUNE_THROW(Dune::InvalidStateException, "Found " << edgeEmbedments.size() << " instead of 2 edge element embedments");
    for (const auto& embedments : edgeEmbedments)
        if (embedments.second.size() != 4)
            DUNE_THROW(Dune::InvalidStateException, "edge element is embedded in " << embedments.second.size() << " facet elements instead of 4");

186
187
    const auto& facetEmbeddings = gridManager.getEmbeddings()->embeddedEntityMap(1);
    const auto& facetEmbedments = gridManager.getEmbeddings()->adjoinedEntityMap(1);
188
189
190
191
192
193
194
195
196
    if (facetEmbeddings.size() != 8) DUNE_THROW(Dune::InvalidStateException, "Found " << facetEmbeddings.size() << " instead of 8 embeddings in facet grid");
    if (facetEmbedments.size() != 32) DUNE_THROW(Dune::InvalidStateException, "Found " << facetEmbedments.size() << " instead of 32 facet element embedments");
    for (const auto& embedments : facetEmbedments)
        if (embedments.second.size() != 2)
            DUNE_THROW(Dune::InvalidStateException, "facet element is embedded in " << embedments.second.size() << " bulk elements instead of 2");
    for (const auto& embeddings : facetEmbeddings)
        if (embeddings.second.size() != 1)
            DUNE_THROW(Dune::InvalidStateException, "facet element has " << embeddings.second.size() << " embedded entities instead of 1");

197
198
    const auto& bulkEmbeddings = gridManager.getEmbeddings()->embeddedEntityMap(0);
    const auto& bulkEmbedments = gridManager.getEmbeddings()->adjoinedEntityMap(0);
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
    if (bulkEmbeddings.size() != 56) DUNE_THROW(Dune::InvalidStateException, "Found " << bulkEmbeddings.size() << " instead of 56 embeddings in bulk grid");
    if (bulkEmbedments.size() != 0) DUNE_THROW(Dune::InvalidStateException, "The grid with highest dimension can't have embedments");

    std::size_t singleEmbeddings = 0;
    std::size_t doubleEmbeddings = 0;
    for (const auto& embeddings : bulkEmbeddings)
    {
        if (embeddings.second.size() != 1 && embeddings.second.size() != 2)
            DUNE_THROW(Dune::InvalidStateException, "bulk element has " << embeddings.second.size() << " embedded entities instead of 1 or 2");
        if (embeddings.second.size() == 1) singleEmbeddings++;
        if (embeddings.second.size() == 2) doubleEmbeddings++;
    }

    if (singleEmbeddings != 48) DUNE_THROW(Dune::InvalidStateException, "Found " << singleEmbeddings << " instead of 48 bulk elements with 1 embedding");
    if (doubleEmbeddings != 8) DUNE_THROW(Dune::InvalidStateException, "Found " << doubleEmbeddings << " instead of 8 bulk elements with 2 embeddings");

    // test access to embeddings for single elements
    singleEmbeddings = 0;
    doubleEmbeddings = 0;
    for (const auto& e : elements(bulkGridView))
    {
220
        if (gridManager.getEmbeddings()->embeddedEntityIndices<0>(e).size() == 1)
221
            singleEmbeddings++;
222
        else if (gridManager.getEmbeddings()->embeddedEntityIndices<0>(e).size() == 2)
223
            doubleEmbeddings++;
224
        else if (gridManager.getEmbeddings()->embeddedEntityIndices<0>(e).size() != 0)
225
226
            DUNE_THROW(Dune::InvalidStateException, "wrong number of embeddings!");

227
        if (gridManager.getEmbeddings()->adjoinedEntityIndices<0>(e).size() != 0)
228
229
230
231
232
233
234
235
236
237
            DUNE_THROW(Dune::InvalidStateException, "bulk grid can't be embedded anywhere!");
    }

    if (singleEmbeddings != 48) DUNE_THROW(Dune::InvalidStateException, "Found " << singleEmbeddings << " instead of 48 bulk elements with 1 embedding");
    if (doubleEmbeddings != 8) DUNE_THROW(Dune::InvalidStateException, "Found " << doubleEmbeddings << " instead of 8 bulk elements with 2 embeddings");

    std::size_t embeddings = 0;
    std::size_t embedments = 0;
    for (const auto& e : elements(facetGridView))
    {
238
        if (gridManager.getEmbeddings()->embeddedEntityIndices<1>(e).size() == 1)
239
            embeddings++;
240
        else if (gridManager.getEmbeddings()->embeddedEntityIndices<1>(e).size() != 0)
241
242
            DUNE_THROW(Dune::InvalidStateException, "wrong number of embeddings!");

243
        if (gridManager.getEmbeddings()->adjoinedEntityIndices<1>(e).size() == 2)
244
            embedments++;
245
        else if (gridManager.getEmbeddings()->adjoinedEntityIndices<1>(e).size() != 0)
246
247
248
249
250
251
252
253
254
255
            DUNE_THROW(Dune::InvalidStateException, "wrong number of embedments!");
    }

    if (embeddings != 8) DUNE_THROW(Dune::InvalidStateException, "Found " << embeddings << " instead of 8 embeddings in facet grid");
    if (embedments != 32) DUNE_THROW(Dune::InvalidStateException, "Found " << embedments << " instead of 32 facet element embedments");

    embeddings = 0;
    embedments = 0;
    for (const auto& e : elements(edgeGridView))
    {
256
        if (gridManager.getEmbeddings()->embeddedEntityIndices<2>(e).size() != 0)
257
258
            DUNE_THROW(Dune::InvalidStateException, "wrong number of embeddings!");

259
        if (gridManager.getEmbeddings()->adjoinedEntityIndices<2>(e).size() == 4)
260
            embedments++;
261
        else if (gridManager.getEmbeddings()->adjoinedEntityIndices<2>(e).size() != 0)
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
            DUNE_THROW(Dune::InvalidStateException, "wrong number of embedments!");
    }

    if (embeddings != 0) DUNE_THROW(Dune::InvalidStateException, "The grid with lowest dimension can't have embedded entities");
    if (embedments!= 2) DUNE_THROW(Dune::InvalidStateException, "Found " << edgeEmbedments.size() << " instead of 2 edge element embedments");

    // write .vtk file for each grid
    using BulkWriter = Dune::VTKWriter<typename BulkGrid::LeafGridView>;
    BulkWriter bulkWriter(bulkGridView); bulkWriter.write("bulkgrid");

    using FacetWriter = Dune::VTKWriter<typename FacetGrid::LeafGridView>;
    FacetWriter facetWriter(facetGridView); facetWriter.write("facetgrid");

    using EdgeWriter = Dune::VTKWriter<typename EdgeGrid::LeafGridView>;
    EdgeWriter edgeWriter(edgeGridView); edgeWriter.write("edgegrid");
}