OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5135 from akva2/griddataoutput_tu

Browse Source 
GridDataOutput: some cleanup
This commit is contained in:
Bård Skaflestad 2024-01-30 14:10:18 +01:00 committed by GitHub
commit 59d781c11b
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
6 changed files with 623 additions and 448 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
CMakeLists_files.cmake
View File

@ -621,8 +621,10 @@ if (Damaris_FOUND AND MPI_FOUND)
list (APPEND PUBLIC_HEADER_FILES ebos/damaris_properties.hh)
list (APPEND PUBLIC_HEADER_FILES ebos/damariswriter.hh)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/DamarisVar.hpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/GridDataOutput.hpp)
list(APPEND MAIN_SOURCE_FILES opm/simulators/utils/DamarisVar.cpp)
list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/GridDataOutput.hpp
opm/simulators/utils/GridDataOutput_impl.hpp)
list(APPEND MAIN_SOURCE_FILES opm/simulators/utils/DamarisVar.cpp
opm/simulators/utils/GridDataOutput.cpp)
endif()
if(HDF5_FOUND)

1
flow/flow_blackoil_alugrid.cpp
View File

@ -35,6 +35,7 @@
#include <ebos/eclgenericwriter_impl.hh>
#include <ebos/ecltransmissibility_impl.hh>
#include <ebos/equil/initstateequil_impl.hh>
#include <opm/simulators/utils/GridDataOutput_impl.hpp>
namespace Opm {
namespace Properties {

1
flow/flow_blackoil_polyhedralgrid.cpp
View File

@ -33,6 +33,7 @@
#include <ebos/ecltransmissibility_impl.hh>
#include <ebos/eclgenericwriter_impl.hh>
#include <ebos/equil/initstateequil_impl.hh>
#include <opm/simulators/utils/GridDataOutput_impl.hpp>
namespace Opm {
namespace Properties {

54
opm/simulators/utils/GridDataOutput.cpp Normal file
View File

@ -0,0 +1,54 @@
/*
Copyright 2023 Inria, Bretagne–Atlantique Research Center
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/utils/GridDataOutput_impl.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/simulators/utils/DamarisVar.hpp>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
#include <dune/fem/gridpart/common/gridpart2gridview.hh>
#include <ebos/femcpgridcompat.hh>
#endif // HAVE_DUNE_FEM
namespace Opm::GridDataOutput {
template<class T> using DV = DamarisOutput::DamarisVar<T>;
#define INSTANCE(part, ...) \
template class SimMeshDataAccessor<__VA_ARGS__, part>; \
template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
writeGridPoints<DV<double>>(DV<double>&, DV<double>&, DV<double>&) const; \
template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
writeConnectivity<DV<int>>(DV<int>&, ConnectivityVertexOrder) const; \
template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
writeOffsetsCells<DV<int>>(DV<int>&) const; \
template long SimMeshDataAccessor<__VA_ARGS__,part>:: \
writeCellTypes<DV<char>>(DV<char>&) const;
INSTANCE(1, Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>)
#if HAVE_DUNE_FEM
INSTANCE(1, Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, (Dune::PartitionIteratorType)4, false> >)
INSTANCE(1, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>)
#endif
} // namespace Opm::GridDataOutput

483
opm/simulators/utils/GridDataOutput.hpp
View File

@ -22,11 +22,10 @@
#include <opm/common/ErrorMacros.hpp>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/grid/io/file/vtk/common.hh>
#include <dune/grid/common/partitionset.hh>
#include <cstddef>
#include <sstream>
#include <iosfwd>
#include <string>
/** @file
@brief Allows model geometry data to be passed to external code - via a copy
@ -129,14 +128,8 @@ public:
* polyhedralCellPresent() to test if polyhedral cells are present and decide
* what they want to do before copying data using the data accessor methods.
*/
explicit SimMeshDataAccessor(const GridView& gridView, Dune::PartitionSet<partitions> dunePartition)
: gridView_(gridView)
, dunePartition_(dunePartition)
{
dimw_ = GridView::dimension; // this is an enum
partition_value_ = dunePartition.value;
countEntities();
}
explicit SimMeshDataAccessor(const GridView& gridView,
Dune::PartitionSet<partitions> dunePartition);
/**
Checks for cells that have polyhedral type within the current partition of
@ -146,16 +139,7 @@ public:
partition is not going to be available for visualisation as this class does
not yet handle polyhedral cells.
*/
bool polyhedralCellPresent()
{
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto corner_geom = cit.geometry();
if (Dune::VTK::geometryType(corner_geom.type()) == Dune::VTK::polyhedron) {
return true;
}
}
return false;
}
bool polyhedralCellPresent() const;
/**
Count the vertices, cells and corners.
@ -164,22 +148,7 @@ public:
do not need to renumber the vertices as all the subsets use references to
the full set.
*/
void countEntities()
{
// We include all the vertices for this ranks partition
const auto& vert_partition = vertices(gridView_, Dune::Partitions::all);
nvertices_ = std::distance(vert_partition.begin(), vert_partition.end());
const auto& cell_partition = elements(gridView_, dunePartition_);
ncells_ = 0;
ncorners_ = 0;
for (const auto& cit : cell_partition) {
auto corner_geom = cit.geometry();
ncorners_ += corner_geom.corners();
++ncells_;
}
}
void countEntities();
/**
Write the positions of vertices - directly to the pointers given in
@ -194,38 +163,7 @@ public:
Returns the number of vertices written
*/
template <typename T>
long writeGridPoints(T* x_inout, T* y_inout, T* z_inout, long max_size = 0)
{
if (max_size < nvertices_) {
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints( T& x_inout, T& "
"y_inout, T& z_inout ) "
+ " Input objects max_size (" + std::to_string(max_size)
+ ") is not sufficient to fit the nvertices_ values (" + std::to_string(nvertices_) + ")");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout[i] = static_cast<T>(xyz_local[0]);
y_inout[i] = static_cast<T>(xyz_local[1]);
z_inout[i] = static_cast<T>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout[i] = static_cast<T>(xyz_local[0]);
y_inout[i] = static_cast<T>(xyz_local[1]);
z_inout[i] = static_cast<T>(0.0);
i++;
}
}
assert(i == nvertices_); // As we are templated on the
// Dune::PartitionSet<partitions>, this cannot change
return i;
}
long writeGridPoints(T* x_inout, T* y_inout, T* z_inout, long max_size = 0) const;
/**
Write the positions of vertices - directly to the pointers given in
@ -241,49 +179,7 @@ public:
Returns the number of vertices written
*/
template <typename VectType>
long writeGridPoints(VectType& x_inout, VectType& y_inout, VectType& z_inout)
{
const std::size_t check_size_x = x_inout.size();
const std::size_t check_size_y = y_inout.size();
const std::size_t check_size_z = z_inout.size();
using VT = decltype(x_inout.data()[0]);
if ((check_size_x < static_cast<std::size_t>(nvertices_)) ||
(check_size_y < static_cast<std::size_t>(nvertices_)) ||
(check_size_z < static_cast<std::size_t>(nvertices_))) {
// assert(check_size >= nvertices_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints( VectType& x_inout, VectType& "
"y_inout, VectType& z_inout ) At least one of the inputs"
+ " object x size " + std::to_string(check_size_x) + " object y size "
+ std::to_string(check_size_y) + " object z size " + std::to_string(check_size_z)
+ " is not sufficient to fit the nvertices_ values( " + std::to_string(nvertices_) + " )");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
z_inout.data()[i] = static_cast<VT>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
double td = 0.0;
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
z_inout.data()[i] = static_cast<VT>(td);
i++;
}
}
assert(i == nvertices_); // As we are templated on the
// Dune::PartitionSet<partitions>, this cannot change
return i;
}
long writeGridPoints(VectType& x_inout, VectType& y_inout, VectType& z_inout) const;
/**
Write the positions of vertices - directly to the pointers given in
@ -296,35 +192,7 @@ public:
Returns the number of vertices written
*/
template <typename T>
long writeGridPoints_AOS(T* xyz_inout, long max_size = 0)
{
if (max_size < nvertices_ * 3) {
assert(max_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_AOS( T* xyz_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
+ std::to_string(nvertices_ * 3) + ")");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i++] = static_cast<T>(xyz_local[0]);
xyz_inout[i++] = static_cast<T>(xyz_local[1]);
xyz_inout[i++] = static_cast<T>(xyz_local[2]);
}
} else if (dimw_ == 2) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i++] = static_cast<T>(xyz_local[0]);
xyz_inout[i++] = static_cast<T>(xyz_local[1]);
xyz_inout[i++] = static_cast<T>(0.0);
}
}
return ((i) / 3);
}
long writeGridPoints_AOS(T* xyz_inout, long max_size = 0) const;
/**
Write the positions of vertices - directly to the pointers given in
@ -336,40 +204,7 @@ public:
Returns the number of vertices written
*/
template <typename VectType>
long writeGridPoints_AOS(VectType& xyz_inout)
{
const std::size_t check_size = xyz_inout.size();
using VT = decltype(xyz_inout.data()[0]);
if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
assert(check_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_AOS( VectType& xyz_inout ) "
+ " Input objects check_size (" + std::to_string(check_size)
+ ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
+ ")");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout.data()[i++] = static_cast<VT>(xyz_local[0]);
xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
xyz_inout[i++] = static_cast<VT>(xyz_local[2]);
}
} else if (dimw_ == 2) {
double td = 0.0;
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i++] = static_cast<VT>(xyz_local[0]);
xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
xyz_inout[i++] = static_cast<VT>(td);
}
}
return ((i) / 3);
}
long writeGridPoints_AOS(VectType& xyz_inout) const;
/**
Write the positions of vertices - directly to the pointers given in
@ -382,40 +217,7 @@ public:
Returns the number of vertices written
*/
template <typename T>
long writeGridPoints_SOA(T* xyz_inout, long max_size = 0)
{
if (max_size < nvertices_ * 3) {
// assert(max_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_SOA( T& xyz_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
+ std::to_string(nvertices_ * 3) + ")");
}
long i = 0;
// Get offsets into structure
T* xyz_inout_y = xyz_inout + nvertices_;
T* xyz_inout_z = xyz_inout + (2 * nvertices_);
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i] = static_cast<T>(xyz_local[0]);
xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
xyz_inout_z[i] = static_cast<T>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i] = static_cast<T>(xyz_local[0]);
xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
xyz_inout_z[i] = static_cast<T>(0.0);
i++;
}
}
return (i);
}
long writeGridPoints_SOA(T* xyz_inout, long max_size = 0) const;
/**
Write the positions of vertices - directly to the pointers given in
@ -427,46 +229,7 @@ public:
Returns the number of vertices written
*/
template <typename VectType>
long writeGridPoints_SOA(VectType& xyz_inout)
{
const std::size_t check_size = xyz_inout.size();
if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
// assert(check_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_SOA( VectType& xyz_inout ) "
+ " Input objects check_size (" + std::to_string(check_size)
+ ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
+ ")");
}
using VT = decltype(xyz_inout.data()[0]);
long i = 0;
// Get offsets into structure
VT* xyz_inout_y = xyz_inout.data() + nvertices_;
VT* xyz_inout_z = xyz_inout.data() + (2 * nvertices_);
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
xyz_inout_z[i] = static_cast<VT>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
double td = 0.0;
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
xyz_inout_z[i] = static_cast<VT>(td);
i++;
}
}
return (i);
}
long writeGridPoints_SOA(VectType& xyz_inout) const;
/**
* Write the connectivity array - directly to the pointer given in parameter 1
@ -480,42 +243,8 @@ public:
Returns the number of corner indices written.
*/
template <typename Integer>
long writeConnectivity(Integer* connectivity_inout, ConnectivityVertexOrder whichOrder, long max_size = 0)
{
if (max_size < ncorners_) {
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeConnectivity( T* connectivity_inout,... ) "
+ " Input max_size value (" + std::to_string(max_size)
+ ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
}
long i = 0;
if (whichOrder == DUNE) {
// DUNE order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
connectivity_inout[i + vx] = vxIdx;
}
i += cell_corners;
}
} else {
// VTK order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
int vtkOrder;
vtkOrder = Dune::VTK::renumber(cit.type(), vx);
connectivity_inout[i + vtkOrder] = vxIdx;
}
i += cell_corners;
}
}
return (i);
}
long writeConnectivity(Integer* connectivity_inout,
ConnectivityVertexOrder whichOrder, long max_size = 0) const;
/**
* Write the connectivity array - directly to a VectType object given in parameter 1
@ -530,46 +259,8 @@ public:
Returns the number of corner indices written.
*/
template <typename VectType>
long writeConnectivity(VectType& connectivity_inout, ConnectivityVertexOrder whichOrder)
{
const std::size_t check_size = connectivity_inout.size();
if (check_size < static_cast<std::size_t>(ncorners_)) {
// assert(check_size >= ncorners_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeConnectivity( VectType& "
"connectivity_inout ) "
+ " Input objects size (" + std::to_string(check_size)
+ ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
}
long i = 0;
if (whichOrder == DUNE) {
// DUNE order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
connectivity_inout.data()[i + vx] = vxIdx;
}
i += cell_corners;
}
} else {
// VTK order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
int vtkOrder;
vtkOrder = Dune::VTK::renumber(cit.type(), vx);
connectivity_inout.data()[i + vtkOrder] = vxIdx;
}
i += cell_corners;
}
}
return (i);
}
long writeConnectivity(VectType& connectivity_inout,
ConnectivityVertexOrder whichOrder) const;
/**
* Write the offsets values - directly to the pointer given in parameter 1
@ -583,24 +274,7 @@ public:
Returns number of offset values written + 1
*/
template <typename Integer>
long writeOffsetsCells(Integer* offsets_inout, long max_size = 0)
{
if (max_size < ncells_) {
// assert(max_size >= ncells_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeOffsetsCells( T* offsets_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
+ std::to_string(ncells_) + ")");
}
long i = 1;
offsets_inout[0] = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
offsets_inout[i] = offsets_inout[i - 1] + cell_corners;
i++;
}
return (i); // This should be 1 greater than ncells_
}
long writeOffsetsCells(Integer* offsets_inout, long max_size = 0) const;
/**
* Write the offsets values - directly to a VectType object given in parameter 1
@ -613,29 +287,7 @@ public:
Returns number of offset values written + 1
*/
template <typename VectType>
long writeOffsetsCells(VectType& offsets_inout)
{
const std::size_t check_size = offsets_inout.size();
if (check_size < static_cast<std::size_t>(ncells_)) {
// assert(check_size >= ncells_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeOffsetsCells( VectType& "
"offsets_inout ) "
+ " Input objects size (" + std::to_string(offsets_inout.size())
+ ") is not sufficient to fit the ncells_ values (" + std::to_string(ncells_) + ")");
}
// using VT = decltype(offsets_inout.data()[0]);
long i = 1;
offsets_inout.data()[0] = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
offsets_inout.data()[i] = offsets_inout.data()[i - 1] + cell_corners;
i++;
}
return (i); // This should be 1 greater than ncells_
}
long writeOffsetsCells(VectType& offsets_inout) const;
/**
* Write the cell types values - directly to the pointer given in parameter 1
@ -648,22 +300,7 @@ public:
Returns number of cells type values written
*/
template <typename Integer>
long writeCellTypes(Integer* types_inout, long max_size = 0)
{
if (max_size < ncells_) {
// assert(max_size >= ncells_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeCellTypes( T* types_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
+ std::to_string(ncells_) + ")");
}
int i = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
Integer vtktype = static_cast<Integer>(Dune::VTK::geometryType(cit.type()));
types_inout[i++] = vtktype;
}
return (i);
}
long writeCellTypes(Integer* types_inout, long max_size = 0) const;
/**
* Write the cell types values - directly to the VectType object given in parameter 1
@ -674,91 +311,45 @@ public:
Returns number of cells type values written
*/
template <typename VectType>
long writeCellTypes(VectType& types_inout)
{
const std::size_t check_size = types_inout.size();
long writeCellTypes(VectType& types_inout) const;
if (check_size < static_cast<std::size_t>(ncells_)) {
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeCellTypes( VectType& types_inout ) " + " Input objects check_size ("
+ std::to_string(check_size) + ") is not sufficient to fit the ncells_ values ("
+ std::to_string(ncells_) + ")");
}
int i = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
int vtktype = static_cast<int>(Dune::VTK::geometryType(cit.type()));
types_inout.data()[i++] = vtktype;
}
return (i);
}
std::string getPartitionTypeString()
{
if (this->dunePartition_ == Dune::Partitions::all)
return (std::string("Dune::Partitions::all"));
if (this->dunePartition_ == Dune::Partitions::interior)
return (std::string("Dune::Partitions::interior"));
if (this->dunePartition_ == Dune::Partitions::interiorBorder)
return (std::string("Dune::Partitions::interiorBorder"));
if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlap)
return (std::string("Dune::Partitions::interiorBorderOverlap"));
if (this->dunePartition_ == Dune::Partitions::front)
return (std::string("Dune::Partitions::front"));
if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlapFront)
return (std::string("Dune::Partitions::InteriorBorderOverlapFront"));
if (this->dunePartition_ == Dune::Partitions::border)
return (std::string("Dune::Partitions::border"));
if (this->dunePartition_ == Dune::Partitions::ghost)
return (std::string("Dune::Partitions::ghost"));
return (std::string("Unknown Dune::PartitionSet<>"));
}
std::string getPartitionTypeString() const;
Dune::PartitionSet<partitions> getPartition(void)
{
return (this->dunePartition_);
return this->dunePartition_;
}
void printGridDetails(std::ostream& outstr)
void printGridDetails(std::ostream& outstr) const;
int getNCells() const
{
outstr << "Dune Partition = " << partition_value_ << ", " << getPartitionTypeString() << std::endl;
outstr << "ncells_: " << getNCells() << std::endl;
outstr << "nvertices_: " << getNVertices() << std::endl;
outstr << "ncorners_: " << getNCorners() << std::endl;
return ncells_;
}
int getNCells()
int getNVertices() const
{
return (ncells_);
return nvertices_;
}
int getNVertices()
int getNCorners() const
{
return (nvertices_);
return ncorners_;
}
int getNCorners()
std::string getError() const
{
return (ncorners_);
}
std::string getError()
{
return error_strm_.str();
return error_;
}
void clearError()
{
error_strm_.str("");
error_.clear();
}
bool hasError()
bool hasError() const
{
if (error_strm_.str().length() > 0)
return true;
else
return false;
return !error_.empty();
}
protected:
@ -783,7 +374,7 @@ protected:
int dimw_; // dimensions of the input grid
private:
std::stringstream error_strm_;
std::string error_;
};
} // namespace Opm::GridDataOutput

526
opm/simulators/utils/GridDataOutput_impl.hpp Normal file
View File

@ -0,0 +1,526 @@
/*
Copyright 2023 Inria, Bretagne–Atlantique Research Center
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/simulators/utils/GridDataOutput.hpp>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/grid/io/file/vtk/common.hh>
#include <opm/common/ErrorMacros.hpp>
#include <cassert>
#include <cstddef>
#include <iterator>
#include <ostream>
namespace Opm::GridDataOutput {
template <class GridView, unsigned int partitions>
SimMeshDataAccessor<GridView,partitions>::
SimMeshDataAccessor(const GridView& gridView, Dune::PartitionSet<partitions> dunePartition)
: gridView_(gridView)
, dunePartition_(dunePartition)
{
dimw_ = GridView::dimension; // this is an enum
partition_value_ = dunePartition.value;
countEntities();
}
template <class GridView, unsigned int partitions>
bool SimMeshDataAccessor<GridView,partitions>::polyhedralCellPresent() const
{
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto corner_geom = cit.geometry();
if (Dune::VTK::geometryType(corner_geom.type()) == Dune::VTK::polyhedron) {
return true;
}
}
return false;
}
template <class GridView, unsigned int partitions>
void SimMeshDataAccessor<GridView,partitions>::countEntities()
{
// We include all the vertices for this ranks partition
const auto& vert_partition = vertices(gridView_, Dune::Partitions::all);
nvertices_ = std::distance(vert_partition.begin(), vert_partition.end());
const auto& cell_partition = elements(gridView_, dunePartition_);
ncells_ = 0;
ncorners_ = 0;
for (const auto& cit : cell_partition) {
auto corner_geom = cit.geometry();
ncorners_ += corner_geom.corners();
++ncells_;
}
}
template <class GridView, unsigned int partitions>
template <typename T>
long SimMeshDataAccessor<GridView,partitions>::
writeGridPoints(T* x_inout, T* y_inout, T* z_inout, long max_size) const
{
if (max_size < nvertices_) {
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints( T& x_inout, T& "
"y_inout, T& z_inout ) "
+ " Input objects max_size (" + std::to_string(max_size)
+ ") is not sufficient to fit the nvertices_ values (" + std::to_string(nvertices_) + ")");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout[i] = static_cast<T>(xyz_local[0]);
y_inout[i] = static_cast<T>(xyz_local[1]);
z_inout[i] = static_cast<T>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout[i] = static_cast<T>(xyz_local[0]);
y_inout[i] = static_cast<T>(xyz_local[1]);
z_inout[i] = static_cast<T>(0.0);
i++;
}
}
assert(i == nvertices_); // As we are templated on the
// Dune::PartitionSet<partitions>, this cannot change
return i;
}
template <class GridView, unsigned int partitions>
template <typename VectType>
long SimMeshDataAccessor<GridView,partitions>::
writeGridPoints(VectType& x_inout, VectType& y_inout, VectType& z_inout) const
{
const std::size_t check_size_x = x_inout.size();
const std::size_t check_size_y = y_inout.size();
const std::size_t check_size_z = z_inout.size();
using VT = decltype(x_inout.data()[0]);
if ((check_size_x < static_cast<std::size_t>(nvertices_)) ||
(check_size_y < static_cast<std::size_t>(nvertices_)) ||
(check_size_z < static_cast<std::size_t>(nvertices_))) {
// assert(check_size >= nvertices_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints( VectType& x_inout, VectType& "
"y_inout, VectType& z_inout ) At least one of the inputs"
+ " object x size " + std::to_string(check_size_x) + " object y size "
+ std::to_string(check_size_y) + " object z size " + std::to_string(check_size_z)
+ " is not sufficient to fit the nvertices_ values( " + std::to_string(nvertices_) + " )");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
z_inout.data()[i] = static_cast<VT>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
double td = 0.0;
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0); // vertices only have one corner
x_inout.data()[i] = static_cast<VT>(xyz_local[0]);
y_inout.data()[i] = static_cast<VT>(xyz_local[1]);
z_inout.data()[i] = static_cast<VT>(td);
i++;
}
}
assert(i == nvertices_); // As we are templated on the
// Dune::PartitionSet<partitions>, this cannot change
return i;
}
template <class GridView, unsigned int partitions>
template <typename T>
long SimMeshDataAccessor<GridView,partitions>::
writeGridPoints_AOS(T* xyz_inout, long max_size) const
{
if (max_size < nvertices_ * 3) {
assert(max_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_AOS( T* xyz_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
+ std::to_string(nvertices_ * 3) + ")");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i++] = static_cast<T>(xyz_local[0]);
xyz_inout[i++] = static_cast<T>(xyz_local[1]);
xyz_inout[i++] = static_cast<T>(xyz_local[2]);
}
} else if (dimw_ == 2) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i++] = static_cast<T>(xyz_local[0]);
xyz_inout[i++] = static_cast<T>(xyz_local[1]);
xyz_inout[i++] = static_cast<T>(0.0);
}
}
return i / 3;
}
template <class GridView, unsigned int partitions>
template <typename VectType>
long SimMeshDataAccessor<GridView,partitions>::
writeGridPoints_AOS(VectType& xyz_inout) const
{
const std::size_t check_size = xyz_inout.size();
using VT = decltype(xyz_inout.data()[0]);
if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
assert(check_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_AOS( VectType& xyz_inout ) "
+ " Input objects check_size (" + std::to_string(check_size)
+ ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
+ ")");
}
long i = 0;
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout.data()[i++] = static_cast<VT>(xyz_local[0]);
xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
xyz_inout[i++] = static_cast<VT>(xyz_local[2]);
}
} else if (dimw_ == 2) {
double td = 0.0;
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i++] = static_cast<VT>(xyz_local[0]);
xyz_inout[i++] = static_cast<VT>(xyz_local[1]);
xyz_inout[i++] = static_cast<VT>(td);
}
}
return i / 3;
}
template <class GridView, unsigned int partitions>
template <typename T>
long SimMeshDataAccessor<GridView,partitions>::
writeGridPoints_SOA(T* xyz_inout, long max_size) const
{
if (max_size < nvertices_ * 3) {
// assert(max_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_SOA( T& xyz_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the nvertices_ * 3 values ("
+ std::to_string(nvertices_ * 3) + ")");
}
long i = 0;
// Get offsets into structure
T* xyz_inout_y = xyz_inout + nvertices_;
T* xyz_inout_z = xyz_inout + (2 * nvertices_);
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i] = static_cast<T>(xyz_local[0]);
xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
xyz_inout_z[i] = static_cast<T>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout[i] = static_cast<T>(xyz_local[0]);
xyz_inout_y[i] = static_cast<T>(xyz_local[1]);
xyz_inout_z[i] = static_cast<T>(0.0);
i++;
}
}
return i;
}
template <class GridView, unsigned int partitions>
template <typename VectType>
long SimMeshDataAccessor<GridView,partitions>::
writeGridPoints_SOA(VectType& xyz_inout) const
{
const std::size_t check_size = xyz_inout.size();
if (check_size < static_cast<std::size_t>(nvertices_ * 3)) {
// assert(check_size >= nvertices_ * 3);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeGridPoints_SOA( VectType& xyz_inout ) "
+ " Input objects check_size (" + std::to_string(check_size)
+ ") is not sufficient to fit the nvertices_ * 3 values (" + std::to_string(nvertices_ * 3)
+ ")");
}
using VT = decltype(xyz_inout.data()[0]);
long i = 0;
// Get offsets into structure
VT* xyz_inout_y = xyz_inout.data() + nvertices_;
VT* xyz_inout_z = xyz_inout.data() + (2 * nvertices_);
if (dimw_ == 3) {
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
xyz_inout_z[i] = static_cast<VT>(xyz_local[2]);
i++;
}
} else if (dimw_ == 2) {
double td = 0.0;
for (const auto& vit : vertices(gridView_, Dune::Partitions::all)) {
auto xyz_local = vit.geometry().corner(0);
xyz_inout.data()[i] = static_cast<VT>(xyz_local[0]);
xyz_inout_y[i] = static_cast<VT>(xyz_local[1]);
xyz_inout_z[i] = static_cast<VT>(td);
i++;
}
}
return i;
}
template <class GridView, unsigned int partitions>
template <typename Integer>
long SimMeshDataAccessor<GridView,partitions>::
writeConnectivity(Integer* connectivity_inout,
ConnectivityVertexOrder whichOrder, long max_size) const
{
if (max_size < ncorners_) {
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeConnectivity( T* connectivity_inout,... ) "
+ " Input max_size value (" + std::to_string(max_size)
+ ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
}
long i = 0;
if (whichOrder == DUNE) {
// DUNE order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
connectivity_inout[i + vx] = vxIdx;
}
i += cell_corners;
}
} else {
// VTK order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
int vtkOrder;
vtkOrder = Dune::VTK::renumber(cit.type(), vx);
connectivity_inout[i + vtkOrder] = vxIdx;
}
i += cell_corners;
}
}
return i;
}
template <class GridView, unsigned int partitions>
template <typename VectType>
long SimMeshDataAccessor<GridView,partitions>::
writeConnectivity(VectType& connectivity_inout,
ConnectivityVertexOrder whichOrder) const
{
const std::size_t check_size = connectivity_inout.size();
if (check_size < static_cast<std::size_t>(ncorners_)) {
// assert(check_size >= ncorners_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeConnectivity( VectType& "
"connectivity_inout ) "
+ " Input objects size (" + std::to_string(check_size)
+ ") is not sufficient to fit the ncorners_ values (" + std::to_string(ncorners_) + ")");
}
long i = 0;
if (whichOrder == DUNE) {
// DUNE order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
connectivity_inout.data()[i + vx] = vxIdx;
}
i += cell_corners;
}
} else {
// VTK order
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
for (auto vx = 0; vx < cell_corners; ++vx) {
const int vxIdx = gridView_.indexSet().subIndex(cit, vx, 3);
int vtkOrder;
vtkOrder = Dune::VTK::renumber(cit.type(), vx);
connectivity_inout.data()[i + vtkOrder] = vxIdx;
}
i += cell_corners;
}
}
return i;
}
template <class GridView, unsigned int partitions>
template <typename Integer>
long SimMeshDataAccessor<GridView,partitions>::
writeOffsetsCells(Integer* offsets_inout, long max_size) const
{
if (max_size < ncells_) {
// assert(max_size >= ncells_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeOffsetsCells( T* offsets_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
+ std::to_string(ncells_) + ")");
}
long i = 1;
offsets_inout[0] = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
offsets_inout[i] = offsets_inout[i - 1] + cell_corners;
i++;
}
return i; // This should be 1 greater than ncells_
}
template <class GridView, unsigned int partitions>
template <typename VectType>
long SimMeshDataAccessor<GridView,partitions>::
writeOffsetsCells(VectType& offsets_inout) const
{
const std::size_t check_size = offsets_inout.size();
if (check_size < static_cast<std::size_t>(ncells_)) {
// assert(check_size >= ncells_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeOffsetsCells( VectType& "
"offsets_inout ) "
+ " Input objects size (" + std::to_string(offsets_inout.size())
+ ") is not sufficient to fit the ncells_ values (" + std::to_string(ncells_) + ")");
}
// using VT = decltype(offsets_inout.data()[0]);
long i = 1;
offsets_inout.data()[0] = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
auto cell_corners = cit.geometry().corners();
offsets_inout.data()[i] = offsets_inout.data()[i - 1] + cell_corners;
i++;
}
return i; // This should be 1 greater than ncells_
}
template <class GridView, unsigned int partitions>
template <typename Integer>
long SimMeshDataAccessor<GridView,partitions>::
writeCellTypes(Integer* types_inout, long max_size) const
{
if (max_size < ncells_) {
// assert(max_size >= ncells_);
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeCellTypes( T* types_inout ) " + " Input objects max_size ("
+ std::to_string(max_size) + ") is not sufficient to fit the ncells_ values ("
+ std::to_string(ncells_) + ")");
}
int i = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
Integer vtktype = static_cast<Integer>(Dune::VTK::geometryType(cit.type()));
types_inout[i++] = vtktype;
}
return i;
}
template <class GridView, unsigned int partitions>
template <typename VectType>
long SimMeshDataAccessor<GridView,partitions>::
writeCellTypes(VectType& types_inout) const
{
const std::size_t check_size = types_inout.size();
if (check_size < static_cast<std::size_t>(ncells_)) {
OPM_THROW(std::runtime_error,
"Opm::GridDataOutput::writeCellTypes( VectType& types_inout ) " + " Input objects check_size ("
+ std::to_string(check_size) + ") is not sufficient to fit the ncells_ values ("
+ std::to_string(ncells_) + ")");
}
int i = 0;
for (const auto& cit : elements(gridView_, dunePartition_)) {
int vtktype = static_cast<int>(Dune::VTK::geometryType(cit.type()));
types_inout.data()[i++] = vtktype;
}
return i;
}
template <class GridView, unsigned int partitions>
std::string SimMeshDataAccessor<GridView,partitions>::
getPartitionTypeString() const
{
if (this->dunePartition_ == Dune::Partitions::all) {
return "Dune::Partitions::all";
}
if (this->dunePartition_ == Dune::Partitions::interior) {
return "Dune::Partitions::interior";
}
if (this->dunePartition_ == Dune::Partitions::interiorBorder) {
return "Dune::Partitions::interiorBorder";
}
if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlap) {
return "Dune::Partitions::interiorBorderOverlap";
}
if (this->dunePartition_ == Dune::Partitions::front) {
return "Dune::Partitions::front";
}
if (this->dunePartition_ == Dune::Partitions::interiorBorderOverlapFront) {
return "Dune::Partitions::InteriorBorderOverlapFront";
}
if (this->dunePartition_ == Dune::Partitions::border) {
return "Dune::Partitions::border";
}
if (this->dunePartition_ == Dune::Partitions::ghost) {
return "Dune::Partitions::ghost";
}
return "Unknown Dune::PartitionSet<>";
}
template <class GridView, unsigned int partitions>
void SimMeshDataAccessor<GridView,partitions>::
printGridDetails(std::ostream& outstr) const
{
outstr << "Dune Partition = " << partition_value_ << ", " << getPartitionTypeString() << std::endl;
outstr << "ncells_: " << getNCells() << std::endl;
outstr << "nvertices_: " << getNVertices() << std::endl;
outstr << "ncorners_: " << getNCorners() << std::endl;
}
} // namespace Opm::GridDataOutput