diff --git a/src/hdf5-writer.hh b/src/hdf5-writer.hh index 85a08e1e3b74618edfc1303b4a720be88b9653b7..0ee16c1425c65444c0b9a9c2000de4016083a3be 100644 --- a/src/hdf5-writer.hh +++ b/src/hdf5-writer.hh @@ -20,19 +20,39 @@ class HDF5Writer { using LocalVector = typename Vector::block_type; public: - HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates, - VertexBasis const &vertexBasis, Patch const &surface, - Patch const &frictionalBoundary, - ConvexPolyhedron<LocalVector> const &weakPatch) - : file_(file), + HDF5Writer(HDF5::Grouplike &file, + const std::vector<Vector>& vertexCoordinates, + const std::vector<const VertexBasis* >& vertexBases, + const std::vector<Patch>& surfaces, + const std::vector<Patch>& frictionalBoundaries, + const std::vector<ConvexPolyhedron<LocalVector>>& weakPatches) + : bodyCount_(vertexCoordinates.size()), + file_(file), iterationWriter_(file_), timeWriter_(file_), #if MY_DIM == 3 - patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch), + patchInfoWriters_(bodyCount_), #endif - surfaceWriter_(file_, vertexCoordinates, surface), - frictionalBoundaryWriter_(file_, vertexCoordinates, - frictionalBoundary) { + surfaceWriters_(bodyCount_), + frictionalBoundaryWriters_(bodyCount_) { + + for (size_t i=0; i<bodyCount_; i++) { +#if MY_DIM == 3 + patchInfoWriters_[i] = new PatchInfoWriter<ProgramState, VertexBasis, GridView>(file_, *vertexBases[i], frictionalBoundaries[i], weakPatches[i], i); +#endif + surfaceWriters_[i] = new SurfaceWriter<ProgramState, GridView>(file_, vertexCoordinates[i], surfaces[i], i); + frictionalBoundaryWriters_[i] = new FrictionalBoundaryWriter<ProgramState, GridView>(file_, vertexCoordinates[i], frictionalBoundaries[i], i); + } + } + + ~HDF5Writer() { + for (size_t i=0; i<bodyCount_; i++) { +#if MY_DIM == 3 + delete patchInfoWriters_[i]; +#endif + delete surfaceWriters_[i]; + delete frictionalBoundaryWriters_[i]; + } } template <class Friction> @@ -40,11 +60,14 @@ class HDF5Writer { // for the friction coefficient Friction &friction) { timeWriter_.write(programState); + + for (size_t i=0; i<bodyCount_; i++) { #if MY_DIM == 3 - patchInfoWriter_.write(programState); + patchInfoWriters_[i]->write(programState); #endif - surfaceWriter_.write(programState); - frictionalBoundaryWriter_.write(programState, friction); + surfaceWriters_[i]->write(programState); + frictionalBoundaryWriters_[i]->write(programState, friction); + } } void reportIterations(ProgramState const &programState, @@ -53,14 +76,16 @@ class HDF5Writer { } private: + const size_t bodyCount_; + HDF5::Grouplike &file_; IterationWriter iterationWriter_; TimeWriter<ProgramState> timeWriter_; #if MY_DIM == 3 - PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_; + std::vector<PatchInfoWriter<ProgramState, VertexBasis, GridView>* > patchInfoWriters_; #endif - SurfaceWriter<ProgramState, GridView> surfaceWriter_; - FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_; + std::vector<SurfaceWriter<ProgramState, GridView>* > surfaceWriters_; + std::vector<FrictionalBoundaryWriter<ProgramState, GridView>* > frictionalBoundaryWriters_; }; #endif diff --git a/src/hdf5/frictionalboundary-writer.cc b/src/hdf5/frictionalboundary-writer.cc index 08bd86da79df042ff97746887ff157acdc0efb42..968b404f073b7bdb379d784f1cf736a9c8e1f776 100644 --- a/src/hdf5/frictionalboundary-writer.cc +++ b/src/hdf5/frictionalboundary-writer.cc @@ -10,8 +10,9 @@ template <class ProgramState, class GridView> FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter( HDF5::Grouplike &file, Vector const &vertexCoordinates, - Patch const &frictionalBoundary) - : group_(file, "frictionalBoundary"), + Patch const &frictionalBoundary, size_t const id) + : id_(id), + group_(file, "frictionalBoundary"+std::to_string(id_)), frictionalBoundary_(frictionalBoundary), frictionalBoundaryDisplacementWriter_(group_, "displacement", frictionalBoundary.numVertices(), @@ -35,27 +36,24 @@ template <class ProgramState, class GridView> template <class Friction> void FrictionalBoundaryWriter<ProgramState, GridView>::write( ProgramState const &programState, Friction &friction) { - // TODO + auto const frictionalBoundaryDisplacements = - restrictToSurface(programState.u[0], frictionalBoundary_); + restrictToSurface(programState.u[id_], frictionalBoundary_); addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep, frictionalBoundaryDisplacements); auto const frictionalBoundaryVelocities = - // TODO - restrictToSurface(programState.v[0], frictionalBoundary_); + restrictToSurface(programState.v[id_], frictionalBoundary_); addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep, frictionalBoundaryVelocities); auto const frictionalBoundaryStates = - // TODO - restrictToSurface(programState.alpha[0], frictionalBoundary_); + restrictToSurface(programState.alpha[id_], frictionalBoundary_); addEntry(frictionalBoundaryStateWriter_, programState.timeStep, frictionalBoundaryStates); - // TODO - friction.updateAlpha(programState.alpha[0]); - auto const c = friction.coefficientOfFriction(programState.v[0]); + friction.updateAlpha(programState.alpha[id_]); + auto const c = friction.coefficientOfFriction(programState.v[id_]); auto const frictionalBoundaryCoefficient = restrictToSurface(c, frictionalBoundary_); addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep, diff --git a/src/hdf5/frictionalboundary-writer.hh b/src/hdf5/frictionalboundary-writer.hh index 705b7fb21baf88873032eefe793e1290df750515..4372cdad928bb60643b8996d7571dbc7d14f6cc7 100644 --- a/src/hdf5/frictionalboundary-writer.hh +++ b/src/hdf5/frictionalboundary-writer.hh @@ -12,12 +12,14 @@ template <class ProgramState, class GridView> class FrictionalBoundaryWriter { public: FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates, - Patch const &frictionalBoundary); + Patch const &frictionalBoundary, size_t const id); template <class Friction> void write(ProgramState const &programState, Friction &friction); private: + size_t const id_; + HDF5::Group group_; Patch const &frictionalBoundary_; diff --git a/src/hdf5/frictionalboundary-writer_tmpl.cc b/src/hdf5/frictionalboundary-writer_tmpl.cc index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..3f8251ae4db70b1f9244fa28de7dc0f078cea20d 100644 --- a/src/hdf5/frictionalboundary-writer_tmpl.cc +++ b/src/hdf5/frictionalboundary-writer_tmpl.cc @@ -0,0 +1,13 @@ +#include "../explicitvectors.hh" +#include "../explicitgrid.hh" + +#include "../program_state.hh" + +using MyProgramState = ProgramState<Vector, ScalarVector>; +using MyFriction = GlobalFriction<Matrix, Vector>; + +template class FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>; +template void FrictionalBoundaryWriter<MyProgramState, DeformedGrid::LeafGridView>::write( + MyProgramState const &programState, MyFriction &friction); + + diff --git a/src/hdf5/patchinfo-writer.cc b/src/hdf5/patchinfo-writer.cc index fc723fc7e01e0b882ce7b7db1d49d8bb063160fe..82b394ea14a42f6902e95c1051682c5e9d4b3bd0 100644 --- a/src/hdf5/patchinfo-writer.cc +++ b/src/hdf5/patchinfo-writer.cc @@ -68,8 +68,10 @@ template <class ProgramState, class VertexBasis, class GridView> PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter( HDF5::Grouplike &file, VertexBasis const &vertexBasis, Patch const &frictionalBoundary, - ConvexPolyhedron<LocalVector> const &weakPatch) - : group_(file, "weakPatchGrid"), + ConvexPolyhedron<LocalVector> const &weakPatch, + size_t const id) + : id_(id), + group_(file, "weakPatchGrid"+std::to_string(id_)), vertexBasis_(vertexBasis), gridEvaluator_(weakPatch, frictionalBoundary.gridView()), weakPatchGridVelocityWriter_( @@ -88,7 +90,7 @@ template <class ProgramState, class VertexBasis, class GridView> void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write( ProgramState const &programState) { // TODO - BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[0]); + BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[id_]); auto const gridVelocity = gridEvaluator_.evaluate(velocity); addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity); } diff --git a/src/hdf5/patchinfo-writer.hh b/src/hdf5/patchinfo-writer.hh index 82dea3646910d28bd401c2eb5d35ae7465d2e648..ef98981c9743d0265be0aeafaaa6bcbf25e12e56 100644 --- a/src/hdf5/patchinfo-writer.hh +++ b/src/hdf5/patchinfo-writer.hh @@ -36,11 +36,14 @@ class PatchInfoWriter { public: PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis, Patch const &frictionalBoundary, - ConvexPolyhedron<LocalVector> const &weakPatch); + ConvexPolyhedron<LocalVector> const &weakPatch, + size_t const id); void write(ProgramState const &programState); private: + const size_t id_; + HDF5::Group group_; VertexBasis const &vertexBasis_; diff --git a/src/hdf5/surface-writer.cc b/src/hdf5/surface-writer.cc index 207522bdb9cfbdd6158c56d473a90ce4154901e5..651bb7c264a2450ffc981fd958ecf5f310bb4983 100644 --- a/src/hdf5/surface-writer.cc +++ b/src/hdf5/surface-writer.cc @@ -7,8 +7,9 @@ template <class ProgramState, class GridView> SurfaceWriter<ProgramState, GridView>::SurfaceWriter( - HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface) - : group_(file, "surface"), + HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface, size_t const id) + : id_(id), + group_(file, "surface"+std::to_string(id_)), surface_(surface), surfaceDisplacementWriter_(group_, "displacement", surface.numVertices(), Vector::block_type::dimension), @@ -24,13 +25,12 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter( template <class ProgramState, class GridView> void SurfaceWriter<ProgramState, GridView>::write( ProgramState const &programState) { - // TODO - auto const surfaceDisplacements = restrictToSurface(programState.u[0], surface_); + + auto const surfaceDisplacements = restrictToSurface(programState.u[id_], surface_); addEntry(surfaceDisplacementWriter_, programState.timeStep, surfaceDisplacements); - // TODO - auto const surfaceVelocities = restrictToSurface(programState.v[0], surface_); + auto const surfaceVelocities = restrictToSurface(programState.v[id_], surface_); addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities); } diff --git a/src/hdf5/surface-writer.hh b/src/hdf5/surface-writer.hh index e0ed4263edb289b8d91ca81d228f428ec7402957..b4b47a781c94a246ffc2edda9216769515567e66 100644 --- a/src/hdf5/surface-writer.hh +++ b/src/hdf5/surface-writer.hh @@ -12,11 +12,13 @@ template <class ProgramState, class GridView> class SurfaceWriter { public: SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates, - Patch const &surface); + Patch const &surface, size_t const id); void write(ProgramState const &programState); private: + size_t const id_; + HDF5::Group group_; Patch const &surface_; diff --git a/src/hdf5/surface-writer_tmpl.cc b/src/hdf5/surface-writer_tmpl.cc index 516b59a775d7aa4fc049619cf0006c96150460c1..5f498b5fc05b1714153cdd36202d83b3a94890a1 100644 --- a/src/hdf5/surface-writer_tmpl.cc +++ b/src/hdf5/surface-writer_tmpl.cc @@ -5,4 +5,4 @@ using MyProgramState = ProgramState<Vector, ScalarVector>; -template class SurfaceWriter<MyProgramState, LeafGridView>; +template class SurfaceWriter<MyProgramState, DeformedGrid::LeafGridView>; diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc index e0b3a60fe5a8e5a1f9966f0817f08dd79d3a3609..a6d496ec58f17e4eac506f96767d51aba12a452c 100644 --- a/src/multi-body-problem.cc +++ b/src/multi-body-problem.cc @@ -432,9 +432,16 @@ int main(int argc, char *argv[]) { } - + using MyVertexBasis = typename Assembler::VertexBasis; + using MyCellBasis = typename Assembler::CellBasis; std::vector<Vector> vertexCoordinates(leafVertexCounts.size()); + std::vector<const MyVertexBasis* > vertexBases(leafVertexCounts.size()); + std::vector<const MyCellBasis* > cellBases(leafVertexCounts.size()); + for (size_t i=0; i<vertexCoordinates.size(); i++) { + vertexBases[i] = &(assemblers[i]->vertexBasis); + cellBases[i] = &(assemblers[i]->cellBasis); + auto& vertexCoords = vertexCoordinates[i]; vertexCoords.resize(leafVertexCounts[i]); @@ -444,16 +451,15 @@ int main(int argc, char *argv[]) { vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry()); } - using MyVertexBasis = typename Assembler::VertexBasis; + auto dataWriter = writeData ? std::make_unique< - HDF5Writer<MyProgramState, MyVertexBasis, GridView>>( - *dataFile, vertexCoordinates, myAssembler.vertexBasis, - surface, frictionalBoundary, weakPatch) + HDF5Writer<MyProgramState, MyVertexBasis, DefLeafGridView>>( + *dataFile, vertexCoordinates, vertexBases, + surfaces, frictionalBoundaries, weakPatches) : nullptr; - MyVTKWriter<MyVertexBasis, typename Assembler::CellBasis> const vtkWriter( - myAssembler.cellBasis, myAssembler.vertexBasis, "obs"); + const MyVTKWriter<MyVertexBasis, MyCellBasis> vtkWriter(cellBases, vertexBases, "body"); IterationRegister iterationCount; /* diff --git a/src/program_state.hh b/src/program_state.hh index 1aea46f754672cb4ad8aff92265e7b1864aeb95d..68c465b472a0cc5bee3d28f9b4b32b27270b8bfc 100644 --- a/src/program_state.hh +++ b/src/program_state.hh @@ -66,6 +66,12 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState { } } + ~ProgramState() { + for (size_t i=0; i<bodyCount_; i++) { + delete bodies[i]; + } + } + size_t size() const { return bodyCount_; } @@ -85,7 +91,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState { const Body<Vector::block_type::dimension>& body) { using LocalVector = typename Vector::block_type; - using LocalMatrix = typename Matrix::block_type; + //using LocalMatrix = typename Matrix::block_type; auto constexpr dims = LocalVector::dimension; diff --git a/src/vtk.cc b/src/vtk.cc index 082af7714fc3d6e14ccfcf0a1f7ded8c6ce0c79f..df3c103395e22a4a221d2dac87f44c4f1d5ea150 100644 --- a/src/vtk.cc +++ b/src/vtk.cc @@ -9,50 +9,81 @@ #include "vtk.hh" template <class VertexBasis, class CellBasis> -MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter( - CellBasis const &_cellBasis, VertexBasis const &_vertexBasis, - std::string _prefix) - : cellBasis(_cellBasis), vertexBasis(_vertexBasis), prefix(_prefix) {} +BodyVTKWriter<VertexBasis, CellBasis>::BodyVTKWriter( + CellBasis const & cellBasis, VertexBasis const & vertexBasis, + std::string prefix) + : cellBasis_(cellBasis), vertexBasis_(vertexBasis), prefix_(prefix) {} template <class VertexBasis, class CellBasis> template <class Vector, class ScalarVector> -void MyVTKWriter<VertexBasis, CellBasis>::write( +void BodyVTKWriter<VertexBasis, CellBasis>::write( size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha, ScalarVector const &stress) const { Dune::VTKWriter<typename VertexBasis::GridView> writer( - vertexBasis.getGridView()); + vertexBasis_.getGridView()); auto const displacementPointer = std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>( - vertexBasis, u, "displacement"); + vertexBasis_, u, "displacement"); writer.addVertexData(displacementPointer); auto const velocityPointer = std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>( - vertexBasis, v, "velocity"); + vertexBasis_, v, "velocity"); writer.addVertexData(velocityPointer); auto const AlphaPointer = std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>( - vertexBasis, alpha, "Alpha"); + vertexBasis_, alpha, "Alpha"); writer.addVertexData(AlphaPointer); auto const stressPointer = std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>( - cellBasis, stress, "stress"); + cellBasis_, stress, "stress"); writer.addCellData(stressPointer); - std::string const filename = prefix + std::to_string(record); + std::string const filename = prefix_ + "_" + std::to_string(record); writer.write(filename.c_str(), Dune::VTK::appendedraw); } template <class VertexBasis, class CellBasis> -void MyVTKWriter<VertexBasis, CellBasis>::writeGrid() const { +void BodyVTKWriter<VertexBasis, CellBasis>::writeGrid() const { Dune::VTKWriter<typename VertexBasis::GridView> writer( - vertexBasis.getGridView()); + vertexBasis_.getGridView()); - std::string const filename = prefix + "_grid"; + std::string const filename = prefix_ + "_grid"; writer.write(filename.c_str(), Dune::VTK::appendedraw); } +template <class VertexBasis, class CellBasis> +MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases, + std::string prefix): bodyVTKWriters_(vertexBases.size()), prefix_(prefix) { + + for (size_t i=0; i<bodyVTKWriters_.size(); i++) { + bodyVTKWriters_[i] = new BodyVTKWriter<VertexBasis, CellBasis>(*cellBases[i], *vertexBases[i], prefix_+std::to_string(i)); + } +} + +template <class VertexBasis, class CellBasis> +MyVTKWriter<VertexBasis, CellBasis>::~MyVTKWriter() { + for (size_t i=0; i<bodyVTKWriters_.size(); i++) { + delete bodyVTKWriters_[i]; + } +} + +template <class VertexBasis, class CellBasis> +template <class Vector, class ScalarVector> +void MyVTKWriter<VertexBasis, CellBasis>::write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v, + const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const { + + for (size_t i=0; i<bodyVTKWriters_.size(); i++) { + bodyVTKWriters_[i]->write(record, u[i], v[i], alpha[i], stress[i]); + } +} +template <class VertexBasis, class CellBasis> +void MyVTKWriter<VertexBasis, CellBasis>::writeGrids() const { + for (size_t i=0; i<bodyVTKWriters_.size(); i++) { + bodyVTKWriters_[i]->writeGrid(); + } +} #include "vtk_tmpl.cc" diff --git a/src/vtk.hh b/src/vtk.hh index 733c0d8f25eaa9452661ac783de34ef491d1d0dd..ee21f8993b324aba1f83e27743b167b5ed7395e2 100644 --- a/src/vtk.hh +++ b/src/vtk.hh @@ -3,13 +3,14 @@ #include <string> -template <class VertexBasis, class CellBasis> class MyVTKWriter { - CellBasis const &cellBasis; - VertexBasis const &vertexBasis; - std::string const prefix; +template <class VertexBasis, class CellBasis> class BodyVTKWriter { +private: + CellBasis const &cellBasis_; + VertexBasis const &vertexBasis_; + std::string const prefix_; public: - MyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis, + BodyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis, std::string prefix); template <class Vector, class ScalarVector> @@ -18,4 +19,22 @@ template <class VertexBasis, class CellBasis> class MyVTKWriter { void writeGrid() const; }; + +template <class VertexBasis, class CellBasis> class MyVTKWriter { +private: + std::vector<BodyVTKWriter<VertexBasis, CellBasis>* > bodyVTKWriters_; + std::string const prefix_; + +public: + MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases, + std::string prefix); + + ~MyVTKWriter(); + + template <class Vector, class ScalarVector> + void write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v, + const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const; + + void writeGrids() const; +}; #endif diff --git a/src/vtk_tmpl.cc b/src/vtk_tmpl.cc index ecfecf9a482599f8afb191e007e77a460e4bc59b..600e7560688860256002bd92d08d320b94bc4344 100644 --- a/src/vtk_tmpl.cc +++ b/src/vtk_tmpl.cc @@ -8,11 +8,11 @@ #include <dune/fufem/functionspacebases/p0basis.hh> #include <dune/fufem/functionspacebases/p1nodalbasis.hh> -using MyP0Basis = P0Basis<LeafGridView, double>; -using P1Basis = P1NodalBasis<LeafGridView, double>; +using MyP0Basis = P0Basis<DeformedGrid::LeafGridView, double>; +using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>; template class MyVTKWriter<P1Basis, MyP0Basis>; template void MyVTKWriter<P1Basis, MyP0Basis>::write<Vector, ScalarVector>( - size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha, - ScalarVector const &stress) const; + size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v, + const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;