Commit a78e605e authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/internal-dirichlet-constraints' into 'master'

[disc] Enable optional internal Dirichlet constraints

Closes #705

See merge request !1664
parents 06897d0e 54d8ad20
......@@ -149,7 +149,7 @@ public:
}
};
this->asImp_().evalDirichletBoundaries(applyDirichlet);
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -173,7 +173,7 @@ public:
jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
};
this->asImp_().evalDirichletBoundaries(applyDirichlet);
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -195,7 +195,17 @@ public:
res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
};
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
//! Enforce Dirichlet constraints
template<typename ApplyFunction>
void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
{
// enforce Dirichlet boundary conditions
this->asImp_().evalDirichletBoundaries(applyDirichlet);
// take care of internal Dirichlet constraints (if enabled)
this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
}
/*!
......
......@@ -86,6 +86,20 @@ public:
{
res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
}
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
const auto eqIdx,
const auto pvIdx)
{
res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
auto& row = jac[scvI.dofIndex()];
for (auto col = row.begin(); col != row.end(); ++col)
row[col.index()][eqIdx] = 0.0;
row[scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
};
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -96,6 +110,19 @@ public:
{
this->asImp_().bindLocalViews();
this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
const auto eqIdx,
const auto pvIdx)
{
auto& row = jac[scvI.dofIndex()];
for (auto col = row.begin(); col != row.end(); ++col)
row[col.index()][eqIdx] = 0.0;
row[scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
};
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -106,6 +133,25 @@ public:
this->asImp_().bindLocalViews();
const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(this->element());
res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
const auto eqIdx,
const auto pvIdx)
{
res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
};
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
//! Enforce Dirichlet constraints
template<typename ApplyFunction>
void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
{
// currently Dirichlet boundary conditions are weakly enforced for cc-schemes
// so here, we only take care of internal Dirichlet constraints (if enabled)
this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
}
};
......
......@@ -61,6 +61,7 @@ class FVLocalAssemblerBase
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using Element = typename GridView::template Codim<0>::Entity;
static constexpr auto numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
public:
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
......@@ -210,6 +211,30 @@ public:
}
}
/*!
* \brief Enforces Dirichlet constraints if enabled in the problem
*/
template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<P::enableInternalDirichletConstraints(), int> = 0>
void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
{
// enforce Dirichlet constraints strongly by overwriting partial derivatives with 1 or 0
// and set the residual to (privar - dirichletvalue)
for (const auto& scvI : scvs(this->fvGeometry()))
{
if (this->problem().hasInternalDirichletConstraint(this->element(), scvI))
{
const auto dirichletValues = this->problem().internalDirichlet(this->element(), scvI);
// set the Dirichlet conditions in residual and jacobian
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
applyDirichlet(scvI, dirichletValues, eqIdx, eqIdx);
}
}
}
template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<!P::enableInternalDirichletConstraints(), int> = 0>
void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
{}
//! The problem
const Problem& problem() const
{ return assembler_.problem(); }
......
......@@ -246,6 +246,25 @@ public:
"a dirichlet() or a dirichletAtPos() method.");
}
/*!
* \brief If internal Dirichlet contraints are enabled
* Enables / disables internal (non-boundary) Dirichlet constraints. If this is overloaded
* to return true, the assembler calls problem.hasInternalDirichletConstraint(element, scv).
* This means you have to implement the following member function
*
* bool hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const;
*
* which returns a bool signifying whether the dof associated with the element/scv pair is contraint.
* If true is returned for a dof, the assembler calls problem.internalDiririchlet(element, scv).
* This means you have to additionally implement the following member function
*
* PrimaryVariables internalDiririchlet(const Element& element, const SubControlVolume& scv) const;
*
* which returns the enforced Dirichlet values the dof associated with the element/scv pair.
*/
static constexpr bool enableInternalDirichletConstraints()
{ return false; }
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
......
......@@ -152,7 +152,7 @@ public:
};
// incorporate Dirichlet BCs
this->asImp_().evalDirichletBoundaries(applyDirichlet);
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -196,7 +196,7 @@ public:
res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
};
this->asImp_().evalDirichletBoundaries(applyDirichlet);
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -226,6 +226,16 @@ public:
ElementResidualVector evalLocalSourceResidual(const Element& neighbor) const
{ return this->evalLocalSourceResidual(neighbor, implicit ? this->curElemVolVars() : this->prevElemVolVars()); }
//! Enforce Dirichlet constraints
template<typename ApplyFunction>
void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
{
// enforce Dirichlet boundary conditions
this->asImp_().evalDirichletBoundaries(applyDirichlet);
// take care of internal Dirichlet constraints (if enabled)
this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
}
/*!
* \brief Incorporate Dirichlet boundary conditions
* \param applyDirichlet Lambda function for the BC incorporation on an scv
......
......@@ -127,6 +127,20 @@ public:
if (i != id)
this->assembleJacobianCoupling(i, jacRow, res[globalI], gridVariables);
});
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
const auto eqIdx,
const auto pvIdx)
{
res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
auto& row = jacRow[domainId][scvI.dofIndex()];
for (auto col = row.begin(); col != row.end(); ++col)
row[col.index()][eqIdx] = 0.0;
row[scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
};
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -158,6 +172,16 @@ public:
this->asImp_().bindLocalViews();
const auto globalI = this->fvGeometry().fvGridGeometry().elementMapper().index(this->element());
res[globalI] = this->evalLocalResidual()[0]; // forward to the internal implementation
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
const auto eqIdx,
const auto pvIdx)
{
res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
};
this->asImp_().enforceDirichletConstraints(applyDirichlet);
}
/*!
......@@ -243,6 +267,15 @@ public:
const Problem& problem() const
{ return this->assembler().problem(domainId); }
//! Enforce Dirichlet constraints
template<typename ApplyFunction>
void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
{
// currently Dirichlet boundary conditions are weakly enforced for cc-schemes
// so here, we only take care of internal Dirichlet constraints (if enabled)
this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
}
//! return reference to the coupling manager
CouplingManager& couplingManager()
{ return couplingManager_; }
......
add_subdirectory(convergence)
add_subdirectory(pointsources)
add_subdirectory(incompressible)
add_subdirectory(internaldirichlet)
add_subdirectory(compressible)
add_subdirectory(periodicbc)
add_subdirectory(fracture2d3d)
......
......@@ -55,6 +55,7 @@
#include <dumux/assembly/fvassembler.hh>
#include "problem.hh"
#include "../internaldirichlet/problem.hh"
//! Function to write out the scv-wise velocities (overload for mpfa)
template<class FVGridGeometry, class GridVariables, class Sol,
......
......@@ -53,7 +53,6 @@ namespace Dumux {
template<class TypeTag> class OnePTestProblem;
namespace Properties {
// create the type tag nodes
// Create new type tags
namespace TTag {
struct OnePIncompressible { using InheritsFrom = std::tuple<OneP>; };
......@@ -111,10 +110,7 @@ struct Scalar<TypeTag, TTag::OnePIncompressibleTpfaQuad> { using type = Quad; };
/*!
* \ingroup OnePTests
* \brief Test problem for the incompressible one-phase model:
*
* Can be run as <tt>./test_box1pfv</tt> or
* <tt>./test_cc1pfv</tt>
* \brief Test problem for the incompressible one-phase model
*/
template<class TypeTag>
class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
......
dune_add_test(NAME test_1p_internaldirichlet_tpfa
SOURCES ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/main.cc
LABELS porousmediumflow 1p
COMPILE_DEFINITIONS TYPETAG=OnePInternalDirichletTpfa
COMPILE_DEFINITIONS NUMDIFFMETHOD=DiffMethod::analytic
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_1p_internaldirichlet_tpfa-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_tpfa-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_tpfa ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/params.input -Problem.Name test_1p_internaldirichlet_tpfa -Problem.EnableGravity false")
dune_add_test(NAME test_1p_internaldirichlet_box
SOURCES ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/main.cc
LABELS porousmediumflow 1p
COMPILE_DEFINITIONS TYPETAG=OnePInternalDirichletBox
COMPILE_DEFINITIONS NUMDIFFMETHOD=DiffMethod::analytic
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_1p_internaldirichlet_box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_box-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_internaldirichlet_box ${CMAKE_SOURCE_DIR}/test/porousmediumflow/1p/implicit/incompressible/params.input -Problem.Name test_1p_internaldirichlet_box -Problem.EnableGravity false")
// -*- 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 *
* the Free Software Foundation, either version 3 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
* \ingroup OnePTests
* \brief A test for internal Dirichlet constraints
*/
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_INTERNAL_DIRICHLET_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_INTERNAL_DIRICHLET_HH
#include <test/porousmediumflow/1p/implicit/incompressible/problem.hh>
namespace Dumux {
// forward declarations
template<class TypeTag> class OnePTestProblemInternalDirichlet;
namespace Properties {
// Create new type tags
namespace TTag {
struct OnePInternalDirichlet {};
struct OnePInternalDirichletTpfa { using InheritsFrom = std::tuple<OnePInternalDirichlet, OnePIncompressibleTpfa>; };
struct OnePInternalDirichletBox { using InheritsFrom = std::tuple<OnePInternalDirichlet, OnePIncompressibleBox>; };
} // end namespace TTag
// Set the problem type
template<class TypeTag>
struct Problem<TypeTag, TTag::OnePInternalDirichlet>
{ using type = OnePTestProblemInternalDirichlet<TypeTag>; };
} // end namespace Properties
/*!
* \ingroup OnePTests
* \brief A test for internal Dirichlet constraints
*/
template<class TypeTag>
class OnePTestProblemInternalDirichlet : public OnePTestProblem<TypeTag>
{
using ParentType = OnePTestProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NeumannValues = GetPropType<TypeTag, Properties::NumEqVector>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
OnePTestProblemInternalDirichlet(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry)
{}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param globalPos The position of the center of the finite volume
*/
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
values.setAllNeumann();
return values;
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param globalPos The position of the boundary face's integration point in global coordinates
*
* Negative values mean influx.
* E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
NeumannValues neumannAtPos(const GlobalPosition& globalPos) const
{
const auto& gg = this->fvGridGeometry();
if (globalPos[0] < gg.bBoxMin()[0] + eps_)
return NeumannValues(1e3);
else if (globalPos[1] < gg.bBoxMin()[1] + eps_)
return NeumannValues(-1e3);
else
return NeumannValues(0.0);
}
//! Enable internal Dirichlet constraints
static constexpr bool enableInternalDirichletConstraints()
{ return true; }
/*!
* \brief Tag a degree of freedom to carry internal Dirichlet constraints.
* If true is returned for a dof, the equation for this dof is replaced
* by the constraint that its primary variable values must match the
* user-defined values obtained from the function internalDirichlet(),
* which must be defined in the problem.
*
* \param element The finite element
* \param scv The sub-control volume
*/
bool hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const
{
// the pure Neumann problem is only defined up to a constant
// we create a well-posed problem by fixing the pressure at one dof in the middle of the domain
return (scv.dofIndex() == static_cast<std::size_t>(this->fvGridGeometry().numDofs()/2));
}
/*!
* \brief Define the values of internal Dirichlet constraints for a degree of freedom.
* \param element The finite element
* \param scv The sub-control volume
*/
PrimaryVariables internalDirichlet(const Element& element, const SubControlVolume& scv) const
{ return PrimaryVariables(1e5); }
private:
static constexpr Scalar eps_ = 1.5e-7;
};
} // end namespace Dumux
#endif
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="100" NumberOfPoints="121">
<PointData Scalars="p">
<DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
100000 1.24072e+06 -1.04072e+06 100000 2.6305e+06 1.55452e+06 4.13529e+06 3.19796e+06 5.53808e+06 4.69269e+06 6.66643e+06 5.84341e+06
7.45365e+06 6.62732e+06 7.91318e+06 7.04425e+06 8.12384e+06 7.18575e+06 8.19706e+06 7.21704e+06 8.2132e+06 7.22549e+06 -2.4305e+06 -1.35452e+06
100000 2.57578e+06 4.19043e+06 5.39504e+06 6.17581e+06 6.51896e+06 6.28733e+06 6.25484e+06 6.25379e+06 -3.93529e+06 -2.99796e+06 -2.37578e+06
100000 2.57651e+06 4.04084e+06 4.98735e+06 5.43331e+06 5.27427e+06 5.28578e+06 5.28474e+06 -5.33808e+06 -4.49269e+06 -3.99043e+06 -2.37651e+06
100000 1.8383e+06 3.07381e+06 3.95787e+06 4.31876e+06 4.3135e+06 4.31476e+06 -6.46643e+06 -5.64341e+06 -5.19504e+06 -3.84084e+06 -1.6383e+06
100000 1.50469e+06 2.73511e+06 3.33092e+06 3.34248e+06 3.34855e+06 -7.25365e+06 -6.42732e+06 -5.97581e+06 -4.78735e+06 -2.87381e+06 -1.30469e+06
100000 1.54703e+06 2.33338e+06 2.37887e+06 2.39919e+06 -7.71318e+06 -6.84425e+06 -6.31896e+06 -5.23331e+06 -3.75787e+06 -2.53511e+06 -1.34703e+06
100000 1.33109e+06 1.43949e+06 1.50685e+06 -7.92384e+06 -6.98575e+06 -6.08733e+06 -5.07427e+06 -4.11876e+06 -3.13092e+06 -2.13338e+06 -1.13109e+06
100000 623614 760454 -7.99706e+06 -7.01704e+06 -6.05484e+06 -5.08578e+06 -4.1135e+06 -3.14248e+06 -2.17887e+06 -1.23949e+06 -423614
100000 269152 -8.0132e+06 -7.02549e+06 -6.05379e+06 -5.08474e+06 -4.11476e+06 -3.14855e+06 -2.19919e+06 -1.30685e+06 -560454 -69152.5
100000
</DataArray>
</PointData>
<CellData Scalars="process rank">
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0
</DataArray>
</CellData>
<Points>
<DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
0 0 0 0.1 0 0 0 0.1 0 0.1 0.1 0
0.2 0 0 0.2 0.1 0 0.3 0 0 0.3 0.1 0
0.4 0 0 0.4 0.1 0 0.5 0 0 0.5 0.1 0
0.6 0 0 0.6 0.1 0 0.7 0 0 0.7 0.1 0
0.8 0 0 0.8 0.1 0 0.9 0 0 0.9 0.1 0
1 0 0 1 0.1 0 0 0.2 0 0.1 0.2 0
0.2 0.2 0 0.3 0.2 0 0.4 0.2 0 0.5 0.2 0
0.6 0.2 0 0.7 0.2 0 0.8 0.2 0 0.9 0.2 0
1 0.2 0 0 0.3 0 0.1 0.3 0 0.2 0.3 0
0.3 0.3 0 0.4 0.3 0 0.5 0.3 0 0.6 0.3 0
0.7 0.3 0 0.8 0.3 0 0.9 0.3 0 1 0.3 0
0 0.4 0 0.1 0.4 0 0.2 0.4 0 0.3 0.4 0
0.4 0.4 0 0.5 0.4 0 0.6 0.4 0 0.7 0.4 0
0.8 0.4 0 0.9 0.4 0 1 0.4 0 0 0.5 0
0.1 0.5 0 0.2 0.5 0 0.3 0.5 0 0.4 0.5 0
0.5 0.5 0 0.6 0.5 0 0.7 0.5 0 0.8 0.5 0
0.9 0.5 0 1 0.5 0 0 0.6 0 0.1 0.6 0
0.2 0.6 0 0.3 0.6 0 0.4 0.6 0 0.5 0.6 0
0.6 0.6 0 0.7 0.6 0 0.8 0.6 0 0.9 0.6 0
1 0.6 0 0 0.7 0 0.1 0.7 0 0.2 0.7 0
0.3 0.7 0 0.4 0.7 0 0.5 0.7 0 0.6 0.7 0
0.7 0.7 0 0.8 0.7 0 0.9 0.7 0 1 0.7 0
0 0.8 0 0.1 0.8 0 0.2 0.8 0 0.3 0.8 0
0.4 0.8 0 0.5 0.8 0 0.6 0.8 0 0.7 0.8 0
0.8 0.8 0 0.9 0.8 0 1 0.8 0 0 0.9 0
0.1 0.9 0 0.2 0.9 0 0.3 0.9 0 0.4 0.9 0
0.5 0.9 0 0.6 0.9 0 0.7 0.9 0 0.8 0.9 0
0.9 0.9 0 1 0.9 0 0 1 0 0.1 1 0
0.2 1 0 0.3 1 0 0.4 1 0 0.5 1 0
0.6 1 0 0.7 1 0 0.8 1 0 0.9 1 0
1 1 0
</DataArray>
</Points>
<Cells>
<DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
0 1 3 2 1 4 5 3 4 6 7 5
6 8 9 7 8 10 11 9 10 12 13 11
12 14 15 13 14 16 17 15 16 18 19 17
18 20 21 19 2 3 23 22 3 5 24 23
5 7 25 24 7 9 26 25 9 11 27 26
11 13 28 27 13 15 29 28 15 17 30 29
17 19 31 30 19 21 32 31 22 23 34 33
23 24 35 34 24 25 36 35 25 26 37 36
26 27 38 37 27 28 39 38 28 29 40 39
29 30 41 40 30 31 42 41 31 32 43 42
33 34 45 44 34 35 46 45 35 36 47 46
36 37 48 47 37 38 49 48 38 39 50 49
39 40 51 50 40 41 52 51 41 42 53 52
42 43 54 53 44 45 56 55 45 46 57 56
46 47 58 57 47 48 59 58 48 49 60 59
49 50 61 60 50 51 62 61 51 52 63 62
52 53 64 63 53 54 65 64 55 56 67 66
56 57 68 67 57 58 69 68 58 59 70 69
59 60 71 70 60 61 72 71 61 62 73 72
62 63 74 73 63 64 75 74 64 65 76 75
66 67 78 77 67 68 79 78 68 69 80 79
69 70 81 80 70 71 82 81 71 72 83 82
72 73 84 83 73 74 85 84 74 75 86 85
75 76 87 86 77 78 89 88 78 79 90 89
79 80 91 90 80 81 92 91 81 82 93 92
82 83 94 93 83 84 95 94 84 85 96 95
85 86 97 96 86 87 98 97 88 89 100 99
89 90 101 100 90 91 102 101 91 92 103 102
92 93 104 103 93 94 105 104 94 95 106 105
95 96 107 106 96 97 108 107 97 98 109 108
99 100 111 110 100 101 112 111 101 102 113 112
102 103 114 113 103 104 115 114 104 105 116 115
105 106 117 116 106 107 118 117 107 108 119 118
108 109 120 119
</DataArray>
<DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
4 8 12 16 20 24 28 32 36 40 44 48
52 56 60 64 68 72 76 80 84 88 92 96
100 104 108 112 116 120 124 128 132 136 140 144
148 152 156 160 164 168 172 176 180 184 188 192
196 200 204 208 212 216 220 224 228 232 236 240
244 248 252 256 260 264 268 272 276 280 284 288
292 296 300 304 308 312 316 320 324 328 332 336
340 344 348 352 356 360 364 368 372 376 380 384
388 392 396 400
</DataArray>
<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9 9 9 9
9 9 9 9
</DataArray>
</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
<UnstructuredGrid>
<Piece NumberOfCells="100" NumberOfPoints="121">
<CellData Scalars="p">