From 5c83aa9c7de023c00b09fd8867f802a319e0f0b9 Mon Sep 17 00:00:00 2001 From: Tim Spain Date: Mon, 11 Aug 2025 12:09:15 +0200 Subject: [PATCH 1/5] Read stress components from a restart file. --- core/src/include/gridNames.hpp | 3 +++ .../modules/DynamicsModule/BBMDynamics.cpp | 19 ++++++++++++++++++- .../DynamicsModule/include/BBMDynamics.hpp | 2 ++ core/src/modules/include/IDynamics.hpp | 13 +++++++++++++ .../src/include/BrittleCGDynamicsKernel.hpp | 6 ++++++ 5 files changed, 42 insertions(+), 1 deletion(-) diff --git a/core/src/include/gridNames.hpp b/core/src/include/gridNames.hpp index e3313f89a..76b1f423b 100644 --- a/core/src/include/gridNames.hpp +++ b/core/src/include/gridNames.hpp @@ -25,6 +25,9 @@ static const std::string shearName = "shear"; static const std::string divergenceName = "divergence"; static const std::string sigmaIName = "sigmaI"; static const std::string sigmaIIName = "sigmaII"; +static const std::string stress11Name = "stress11"; +static const std::string stress12Name = "stress12"; +static const std::string stress22Name = "stress22"; static const std::string uWindName = "uwind"; static const std::string vWindName = "vwind"; diff --git a/core/src/modules/DynamicsModule/BBMDynamics.cpp b/core/src/modules/DynamicsModule/BBMDynamics.cpp index bffeb7f4a..185ef43c6 100644 --- a/core/src/modules/DynamicsModule/BBMDynamics.cpp +++ b/core/src/modules/DynamicsModule/BBMDynamics.cpp @@ -10,7 +10,11 @@ namespace Nextsim { static const std::vector namedFields = { uName, vName }; -static const std::map> defaultFields = {}; +static const std::map> defaultFields = { + { stress11Name, { ModelArray::Type::DG, 0. } }, + { stress12Name, { ModelArray::Type::DG, 0. } }, + { stress22Name, { ModelArray::Type::DG, 0. } }, +}; // TODO: We should use getName() here, but it isn't static. static const std::string prefix = "BBMDynamics"; // MEVPDynamics::getName(); @@ -162,6 +166,19 @@ void BBMDynamics::update(const TimestepTime& tst) sigmaII = kernel.getDG0Data(sigmaIIName); } +ModelState BBMDynamics::getStatePrognostic() const +{ + ModelState state = { { + { stress11Name, stress11 }, + { stress12Name, stress12 }, + { stress22Name, stress22 }, + }, + getConfiguration() + }; + state.merge(IDynamics::getStatePrognostic()); + return state; +} + void BBMDynamics::prepareAdvection() { kernel.prepareAdvection(); } void BBMDynamics::advectField( diff --git a/core/src/modules/DynamicsModule/include/BBMDynamics.hpp b/core/src/modules/DynamicsModule/include/BBMDynamics.hpp index a72705ea4..05defe5de 100644 --- a/core/src/modules/DynamicsModule/include/BBMDynamics.hpp +++ b/core/src/modules/DynamicsModule/include/BBMDynamics.hpp @@ -21,6 +21,8 @@ class BBMDynamics : public IDynamics, public Configured { public: BBMDynamics(); + ModelState getStatePrognostic() const override; + std::string getName() const override { return "BBMDynamics"; } void prepareAdvection() override; void update(const TimestepTime& tst) override; diff --git a/core/src/modules/include/IDynamics.hpp b/core/src/modules/include/IDynamics.hpp index 2be8c7ffd..faf2e9c19 100644 --- a/core/src/modules/include/IDynamics.hpp +++ b/core/src/modules/include/IDynamics.hpp @@ -27,6 +27,9 @@ class IDynamics : public ModelComponent { , divergence(ModelArray::Type::H) , sigmaI(ModelArray::Type::H) , sigmaII(ModelArray::Type::H) + , stress11(ModelArray::Type::DG) + , stress12(ModelArray::Type::DG) + , stress22(ModelArray::Type::DG) , damage(getStore()) , uwind(getStore()) , vwind(getStore()) @@ -76,6 +79,9 @@ class IDynamics : public ModelComponent { { divergenceName, divergence }, { sigmaIName, sigmaI }, { sigmaIIName, sigmaII }, + { stress11Name, stress11 }, + { stress11Name, stress12 }, + { stress11Name, stress22 }, }, {} }; return state.merge(getStatePrognostic()); @@ -90,6 +96,9 @@ class IDynamics : public ModelComponent { divergence.resize(); sigmaI.resize(); sigmaII.resize(); + stress11.resize(); + stress12.resize(); + stress22.resize(); } virtual void update(const TimestepTime& tst) = 0; @@ -119,6 +128,10 @@ class IDynamics : public ModelComponent { HField divergence; HField sigmaI; HField sigmaII; + // Stress components. Diagnostic for some dynamics, prognostic for brittle and others. + DGField stress11; + DGField stress12; + DGField stress22; // References to the DG0 finite volume data arrays ModelArrayRef damage; diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index eb80724b0..b46f1b2bb 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -153,6 +153,12 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern { if (name == damageName) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); + } else if (name == stress11Name) { + DGModelArray::ma2dg(data, s11); + } else if (name == stress12Name) { + DGModelArray::ma2dg(data, s12); + } else if (name == stress22Name) { + DGModelArray::ma2dg(data, s22); } else { CGDynamicsKernel::setData(name, data); } From 809f59008554faec5b7d2a2eb1b5432686b1e988 Mon Sep 17 00:00:00 2001 From: Tim Spain Date: Mon, 1 Sep 2025 09:34:35 +0200 Subject: [PATCH 2/5] Add exceptions to catch mismatches in the number of components. --- dynamics/src/include/dgVectorHolder.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/dynamics/src/include/dgVectorHolder.hpp b/dynamics/src/include/dgVectorHolder.hpp index 8b87a3160..b28667a93 100644 --- a/dynamics/src/include/dgVectorHolder.hpp +++ b/dynamics/src/include/dgVectorHolder.hpp @@ -8,6 +8,8 @@ #include "dgVector.hpp" #include "include/ModelArray.hpp" +#include + namespace Nextsim { template class DGVectorHolder { @@ -20,6 +22,10 @@ template class DGVectorHolder { DGVectorHolder(ModelArray& ma) : DGVectorHolder(static_cast(ma)) { + if (ma.nComponents() != DG) { + std::string msg = "DGVectorHolder<"+std::to_string(DG)+">(ModelArray&): incorrect number of components = " + std::to_string(ma.nComponents()); + throw std::length_error(msg); + } } DGVectorHolder(ModelArray::DataType& madt) : DGVectorHolder(reinterpret_cast(madt)) From 6a934343ab7ba32c30e107841ae74560d3c39402 Mon Sep 17 00:00:00 2001 From: Tim Spain Date: Wed, 24 Sep 2025 06:18:55 +0200 Subject: [PATCH 3/5] Move strain & stress and change stress to a DGVectorHolder. --- .../modules/DynamicsModule/BBMDynamics.cpp | 19 ++++++-- core/src/modules/include/IDynamics.hpp | 12 ++--- dynamics/src/CGDynamicsKernel.cpp | 41 +++++++++++++++- .../src/include/BrittleCGDynamicsKernel.hpp | 24 +++++----- dynamics/src/include/CGDynamicsKernel.hpp | 12 ++--- dynamics/src/include/DynamicsKernel.hpp | 48 +++++++++---------- dynamics/src/include/VPCGDynamicsKernel.hpp | 12 ++--- dynamics/src/include/dgVectorHolder.hpp | 3 ++ 8 files changed, 113 insertions(+), 58 deletions(-) diff --git a/core/src/modules/DynamicsModule/BBMDynamics.cpp b/core/src/modules/DynamicsModule/BBMDynamics.cpp index 185ef43c6..3c38c319e 100644 --- a/core/src/modules/DynamicsModule/BBMDynamics.cpp +++ b/core/src/modules/DynamicsModule/BBMDynamics.cpp @@ -11,9 +11,9 @@ namespace Nextsim { static const std::vector namedFields = { uName, vName }; static const std::map> defaultFields = { - { stress11Name, { ModelArray::Type::DG, 0. } }, - { stress12Name, { ModelArray::Type::DG, 0. } }, - { stress22Name, { ModelArray::Type::DG, 0. } }, +// { stress11Name, { ModelArray::Type::DG, 0. } }, +// { stress12Name, { ModelArray::Type::DG, 0. } }, +// { stress22Name, { ModelArray::Type::DG, 0. } }, }; // TODO: We should use getName() here, but it isn't static. @@ -107,6 +107,16 @@ void BBMDynamics::setData(const ModelState::DataMap& ms) uice = ms.at(uName); vice = ms.at(vName); + if (ms.count(stress11Name)) { + stress11 = ms.at(stress11Name); + stress12 = ms.at(stress12Name); + stress22 = ms.at(stress22Name); + } else { + stress11 = 0.; + stress12 = 0.; + stress22 = 0.; + } + // Set the data in the kernel arrays. // Required data for (const auto& fieldName : namedFields) { @@ -134,6 +144,9 @@ void BBMDynamics::setData(const ModelState::DataMap& ms) kernel.setDGArray(ciceName, ciceDG.allComponents()); kernel.setDGArray(hsnowName, hsnowDG.allComponents()); kernel.setDGArray(damageName, damage.allComponents()); + kernel.setDGArray(stress11Name, stress11); + kernel.setDGArray(stress12Name, stress12); + kernel.setDGArray(stress22Name, stress22); } void BBMDynamics::update(const TimestepTime& tst) diff --git a/core/src/modules/include/IDynamics.hpp b/core/src/modules/include/IDynamics.hpp index faf2e9c19..7714db2a1 100644 --- a/core/src/modules/include/IDynamics.hpp +++ b/core/src/modules/include/IDynamics.hpp @@ -27,9 +27,9 @@ class IDynamics : public ModelComponent { , divergence(ModelArray::Type::H) , sigmaI(ModelArray::Type::H) , sigmaII(ModelArray::Type::H) - , stress11(ModelArray::Type::DG) - , stress12(ModelArray::Type::DG) - , stress22(ModelArray::Type::DG) + , stress11(ModelArray::Type::DGSTRESS) + , stress12(ModelArray::Type::DGSTRESS) + , stress22(ModelArray::Type::DGSTRESS) , damage(getStore()) , uwind(getStore()) , vwind(getStore()) @@ -129,9 +129,9 @@ class IDynamics : public ModelComponent { HField sigmaI; HField sigmaII; // Stress components. Diagnostic for some dynamics, prognostic for brittle and others. - DGField stress11; - DGField stress12; - DGField stress22; + DGSField stress11; + DGSField stress12; + DGSField stress22; // References to the DG0 finite volume data arrays ModelArrayRef damage; diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index b9fe73fde..3d3a06890 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -15,6 +15,7 @@ #include "include/ParametricMap.hpp" #include "include/VectorManipulations.hpp" #include "include/cgVector.hpp" +#include "include/DGModelArray.hpp" #include @@ -41,6 +42,10 @@ void CGDynamicsKernel::initialise( u.resize_by_mesh(*smesh); v.resize_by_mesh(*smesh); + e11.resize_by_mesh(*smesh); + e12.resize_by_mesh(*smesh); + e22.resize_by_mesh(*smesh); + cgH.resize_by_mesh(*smesh); cgA.resize_by_mesh(*smesh); @@ -81,6 +86,20 @@ void CGDynamicsKernel::setData(const std::string& name, const Model } } +template +void CGDynamicsKernel::setDGArray(const std::string& name, ModelArray::DataType& dgData) +{ + if (name == stress11Name) { + s11 = DGVectorHolder(dgData); + } else if (name == stress12Name) { + s12 = DGVectorHolder(dgData); + } else if (name == stress22Name) { + s22 = DGVectorHolder(dgData); + } else { + DynamicsKernel::setDGArray(name, dgData); + } +} + template ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) const { @@ -104,7 +123,27 @@ ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) co DGVector vtmp(*smesh); Nextsim::Interpolations::CG2DG(*smesh, vtmp, vIceOceanStress); return DGModelArray::dg2ma(vtmp, data); - } else { + } else if (name == shearName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma(Tools::Shear(*smesh, e11, e12, e22), data); + } else if (name == divergenceName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, e11, e12, e22), data); + } else if (name == sigmaIName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, + static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22) + ), data); + } else if (name == sigmaIIName) { + ModelArray data(ModelArray::Type::H); + return DGModelArray::dg2ma(Tools::TensorInvII(*smesh, + static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22) + ), data); + } else { return DynamicsKernel::getDG0Data(name); } } diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index b46f1b2bb..7aaca18d9 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -20,12 +20,6 @@ namespace Nextsim { // The brittle momentum solver for CG velocity fields template class BrittleCGDynamicsKernel : public CGDynamicsKernel { protected: - using DynamicsKernel::s11; - using DynamicsKernel::s12; - using DynamicsKernel::s22; - using DynamicsKernel::e11; - using DynamicsKernel::e12; - using DynamicsKernel::e22; using DynamicsKernel::hice; using DynamicsKernel::cice; using DynamicsKernel::smesh; @@ -37,6 +31,12 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern using CGDynamicsKernel::u; using CGDynamicsKernel::v; + using CGDynamicsKernel::s11; + using CGDynamicsKernel::s12; + using CGDynamicsKernel::s22; + using CGDynamicsKernel::e11; + using CGDynamicsKernel::e12; + using CGDynamicsKernel::e22; using CGDynamicsKernel::xGradSeaSurfaceHeight; using CGDynamicsKernel::yGradSeaSurfaceHeight; using CGDynamicsKernel::uAtmos; @@ -153,12 +153,12 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern { if (name == damageName) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); - } else if (name == stress11Name) { - DGModelArray::ma2dg(data, s11); - } else if (name == stress12Name) { - DGModelArray::ma2dg(data, s12); - } else if (name == stress22Name) { - DGModelArray::ma2dg(data, s22); +// } else if (name == stress11Name) { +// DGModelArray::ma2dg(data, s11); +// } else if (name == stress12Name) { +// DGModelArray::ma2dg(data, s12); +// } else if (name == stress22Name) { +// DGModelArray::ma2dg(data, s22); } else { CGDynamicsKernel::setData(name, data); } diff --git a/dynamics/src/include/CGDynamicsKernel.hpp b/dynamics/src/include/CGDynamicsKernel.hpp index 61de16217..bd952c1bf 100644 --- a/dynamics/src/include/CGDynamicsKernel.hpp +++ b/dynamics/src/include/CGDynamicsKernel.hpp @@ -26,14 +26,9 @@ namespace Nextsim { template class CGDynamicsKernel : public DynamicsKernel { protected: - using DynamicsKernel::s11; - using DynamicsKernel::s12; - using DynamicsKernel::s22; - using DynamicsKernel::e11; - using DynamicsKernel::e12; - using DynamicsKernel::e22; using DynamicsKernel::smesh; using DynamicsKernel::dgtransport; + using DynamicsKernel::setDGArray; using typename DynamicsKernel::DataMap; public: @@ -42,6 +37,7 @@ class CGDynamicsKernel : public DynamicsKernel { void initialise(const ModelArray& coords, bool isSpherical, const ModelArray& mask) override; void setData(const std::string& name, const ModelArray& data) override; + void setDGArray(const std::string& name, ModelArray::DataType& dgData) override; ModelArray getDG0Data(const std::string& name) const override; void computeGradientOfSeaSurfaceHeight(const DGVector<1>& seaSurfaceHeight); void prepareIteration(const DataMap& data) override; @@ -67,6 +63,10 @@ class CGDynamicsKernel : public DynamicsKernel { CGVector cgA; CGVector cgH; + //! DGVectorHolders for strain and stress components + DGVector e11, e12, e22; + DGVectorHolder s11, s12, s22; + // CG gradient of the seaSurfaceHeight CGVector xGradSeaSurfaceHeight; CGVector yGradSeaSurfaceHeight; diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index fc228c4e9..21d371c8f 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -57,20 +57,20 @@ template class DynamicsKernel { seaSurfaceHeight.resize_by_mesh(*smesh); - e11.resize_by_mesh(*smesh); - e12.resize_by_mesh(*smesh); - e22.resize_by_mesh(*smesh); - s11.resize_by_mesh(*smesh); - s12.resize_by_mesh(*smesh); - s22.resize_by_mesh(*smesh); +// e11.resize_by_mesh(*smesh); +// e12.resize_by_mesh(*smesh); +// e22.resize_by_mesh(*smesh); +// s11.resize_by_mesh(*smesh); +// s12.resize_by_mesh(*smesh); +// s22.resize_by_mesh(*smesh); // Set initial values to zero. Prognostic fields will be filled from the restart file. - e11.zero(); - e12.zero(); - e22.zero(); - s11.zero(); - s12.zero(); - s22.zero(); +// e11.zero(); +// e12.zero(); +// e22.zero(); +// s11.zero(); +// s12.zero(); +// s22.zero(); } /*! @@ -91,9 +91,16 @@ template class DynamicsKernel { */ virtual void setData(const std::string& name, const ModelArray& data) { - + static const std::set vectorHolderNames = { + hiceName, + ciceName, + hsnowName, + stress11Name, + stress12Name, + stress22Name, + }; // Special cases: hice, cice, (damage, stress) <- not yet implemented - if (name == hiceName || name == ciceName || name == hsnowName) { + if (vectorHolderNames.count(name)) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); } else if (name == sshName) { DGModelArray::ma2dg(data, seaSurfaceHeight); @@ -125,14 +132,6 @@ template class DynamicsKernel { if (name == hiceName || name == ciceName || name == hsnowName) { throw std::runtime_error( std::string("DynamicsKernel::getDG0Data: Use array sharing for ") + name); - } else if (name == shearName) { - return DGModelArray::dg2ma(Tools::Shear(*smesh, e11, e12, e22), data); - } else if (name == divergenceName) { - return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, e11, e12, e22), data); - } else if (name == sigmaIName) { - return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, s11, s12, s22), data); - } else if (name == sigmaIIName) { - return DGModelArray::dg2ma(Tools::TensorInvII(*smesh, s11, s12, s22), data); } else { // Any other named field must exist return ModelArray(ModelArray::component0Type(ModelArray::AdvectionType)); @@ -199,8 +198,9 @@ template class DynamicsKernel { DGVector<1> seaSurfaceHeight; //! Vectors storing strain and stress components - DGVector e11, e12, e22; - DGVector s11, s12, s22; +// DGVector e11, e12, e22; +// DGVector /*s11,*/ s12, s22; +// DGVectorHolder s11, s12, s22; size_t stepNumber = 0; diff --git a/dynamics/src/include/VPCGDynamicsKernel.hpp b/dynamics/src/include/VPCGDynamicsKernel.hpp index 2ac19bcb3..fcdb0e107 100644 --- a/dynamics/src/include/VPCGDynamicsKernel.hpp +++ b/dynamics/src/include/VPCGDynamicsKernel.hpp @@ -19,12 +19,6 @@ namespace Nextsim { // The VP pseudo-timestepping momentum equation solver for CG velocities template class VPCGDynamicsKernel : public CGDynamicsKernel { protected: - using DynamicsKernel::s11; - using DynamicsKernel::s12; - using DynamicsKernel::s22; - using DynamicsKernel::e11; - using DynamicsKernel::e12; - using DynamicsKernel::e22; using DynamicsKernel::hice; using DynamicsKernel::cice; using DynamicsKernel::smesh; @@ -35,6 +29,12 @@ template class VPCGDynamicsKernel : public CGDynamicsKernel::u; using CGDynamicsKernel::v; + using CGDynamicsKernel::s11; + using CGDynamicsKernel::s12; + using CGDynamicsKernel::s22; + using CGDynamicsKernel::e11; + using CGDynamicsKernel::e12; + using CGDynamicsKernel::e22; using CGDynamicsKernel::xGradSeaSurfaceHeight; using CGDynamicsKernel::yGradSeaSurfaceHeight; using CGDynamicsKernel::uAtmos; diff --git a/dynamics/src/include/dgVectorHolder.hpp b/dynamics/src/include/dgVectorHolder.hpp index b28667a93..b41585aeb 100644 --- a/dynamics/src/include/dgVectorHolder.hpp +++ b/dynamics/src/include/dgVectorHolder.hpp @@ -47,6 +47,9 @@ template class DGVectorHolder { void zero() { ref->setZero(); } + auto row(Eigen::Index i) { return static_cast&>(*ref).row(i); } + auto col(Eigen::Index i) { return static_cast&>(*ref).col(i); } + private: EigenDGVector* ref; }; From c76fb3b3a0d1bbe207d77c1ed59211cab55929a5 Mon Sep 17 00:00:00 2001 From: Tim Spain Date: Wed, 24 Sep 2025 09:07:19 +0200 Subject: [PATCH 4/5] clang formatting --- dynamics/src/CGDynamicsKernel.cpp | 151 +++++++++--------- .../src/include/BrittleCGDynamicsKernel.hpp | 43 ++--- dynamics/src/include/DynamicsKernel.hpp | 37 ++--- dynamics/src/include/dgVectorHolder.hpp | 48 ++++-- 4 files changed, 145 insertions(+), 134 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index 3d3a06890..b28dc58d8 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -21,15 +21,15 @@ namespace Nextsim { -template +template CGDynamicsKernel::CGDynamicsKernel(const DynamicsParameters& params) - : baseParams(params) + : baseParams(params) { } -template -void CGDynamicsKernel::initialise( - const ModelArray& coords, bool isSpherical, const ModelArray& mask) +template +void CGDynamicsKernel::initialise(const ModelArray& coords, bool isSpherical, + const ModelArray& mask) { DynamicsKernel::initialise(coords, isSpherical, mask); @@ -68,17 +68,12 @@ void CGDynamicsKernel::initialise( sinOceanAngle = std::sin(radians(baseParams.oceanTurningAngle)); } -template +template void CGDynamicsKernel::setData(const std::string& name, const ModelArray& data) { - const std::map*> targetMap = { - { uName, &u }, - { vName, &v }, - { uWindName, &uAtmos }, - { vWindName, &vAtmos }, - { uOceanName, &uOcean }, - { vOceanName, &vOcean }, - }; + const std::map*> targetMap = { { uName, &u }, { vName, &v }, { + uWindName, &uAtmos }, { vWindName, &vAtmos }, { uOceanName, &uOcean }, { vOceanName, + &vOcean }, }; if (targetMap.count(name)) { ma2cg(data, *targetMap.at(name)); } else { @@ -86,8 +81,9 @@ void CGDynamicsKernel::setData(const std::string& name, const Model } } -template -void CGDynamicsKernel::setDGArray(const std::string& name, ModelArray::DataType& dgData) +template +void CGDynamicsKernel::setDGArray(const std::string& name, + ModelArray::DataType& dgData) { if (name == stress11Name) { s11 = DGVectorHolder(dgData); @@ -100,7 +96,7 @@ void CGDynamicsKernel::setDGArray(const std::string& name, ModelArr } } -template +template ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) const { if (name == uName) { @@ -131,31 +127,29 @@ ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) co return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, e11, e12, e22), data); } else if (name == sigmaIName) { ModelArray data(ModelArray::Type::H); - return DGModelArray::dg2ma(Tools::TensorInvI(*smesh, - static_cast&>(s11), - static_cast&>(s12), - static_cast&>(s22) - ), data); + return DGModelArray::dg2ma( + Tools::TensorInvI(*smesh, static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22)), data); } else if (name == sigmaIIName) { ModelArray data(ModelArray::Type::H); - return DGModelArray::dg2ma(Tools::TensorInvII(*smesh, - static_cast&>(s11), - static_cast&>(s12), - static_cast&>(s22) - ), data); - } else { + return DGModelArray::dg2ma( + Tools::TensorInvII(*smesh, static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22)), data); + } else { return DynamicsKernel::getDG0Data(name); } } -template void CGDynamicsKernel::prepareAdvection() +template void CGDynamicsKernel::prepareAdvection() { dgtransport->prepareAdvection(u, v); } -template +template void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( - const DGVector<1>& seaSurfaceHeight) + const DGVector<1>& seaSurfaceHeight) { // First transform to CG1 Vector and do all computations in CG1 CGVector<1> cgSeasurfaceHeight(*smesh); @@ -177,8 +171,9 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( for (size_t cx = 0; cx < smesh->nx; ++cx, ++eid, ++cg1id) { // get local CG nodes Eigen::Vector loc_cgSSH = { cgSeasurfaceHeight(cg1id), - cgSeasurfaceHeight(cg1id + 1), cgSeasurfaceHeight(cg1id + smesh->nx + 1), - cgSeasurfaceHeight(cg1id + smesh->nx + 1 + 1) }; + cgSeasurfaceHeight(cg1id + 1), cgSeasurfaceHeight( + cg1id + smesh->nx + 1), cgSeasurfaceHeight( + cg1id + smesh->nx + 1 + 1) }; // compute grad Eigen::Vector tx = pmap->dX_SSH[eid] * loc_cgSSH; @@ -209,14 +204,14 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( size_t cg1row = smesh->nx + 1; size_t topleft = smesh->ny * cg1row; for (size_t i = 1; i < smesh->nx; ++i) // bottom / top - { + { uGrad(i) = uGrad(i + cg1row); vGrad(i) = vGrad(i + cg1row); uGrad(topleft + i) = uGrad(topleft + i - cg1row); vGrad(topleft + i) = vGrad(topleft + i - cg1row); } for (size_t i = 1; i < smesh->ny; ++i) // left / right - { + { uGrad(i * cg1row) = uGrad(i * cg1row + 1); vGrad(i * cg1row) = vGrad(i * cg1row + 1); @@ -251,7 +246,7 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( // along lines #pragma omp parallel for for (size_t iy = 0; iy <= smesh->ny; ++iy) // horizontal - { + { size_t icg1 = (smesh->nx + 1) * iy; size_t icg2 = (2 * smesh->nx + 1) * 2 * iy + 1; for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) { @@ -261,7 +256,7 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( } #pragma omp parallel for for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical - { + { size_t icg1 = (smesh->nx + 1) * iy; size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1); for (size_t ix = 0; ix <= smesh->nx; ++ix, ++icg1, icg2 += 2) { @@ -273,22 +268,22 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( // midpoints #pragma omp parallel for for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical - { + { size_t icg1 = (smesh->nx + 1) * iy; size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1) + 1; for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) { xGradSeaSurfaceHeight(icg2) = 0.25 - * (uGrad(icg1) + uGrad(icg1 + 1) + uGrad(icg1 + cg1row) - + uGrad(icg1 + cg1row + 1)); + * (uGrad(icg1) + uGrad(icg1 + 1) + uGrad(icg1 + cg1row) + + uGrad(icg1 + cg1row + 1)); yGradSeaSurfaceHeight(icg2) = 0.25 - * (vGrad(icg1) + vGrad(icg1 + 1) + vGrad(icg1 + cg1row) - + vGrad(icg1 + cg1row + 1)); + * (vGrad(icg1) + vGrad(icg1 + 1) + vGrad(icg1 + cg1row) + + vGrad(icg1 + cg1row + 1)); } } } } -template void CGDynamicsKernel::prepareIteration(const DataMap& data) +template void CGDynamicsKernel::prepareIteration(const DataMap& data) { // interpolate ice height and concentration to local cg Variables Interpolations::DG2CG(*smesh, cgH, data.at(hiceName)); @@ -310,11 +305,11 @@ template void CGDynamicsKernel::prepareIteration( cgH = cgH.cwiseMax(0.05); } -template -Eigen::Matrix cgLocal( - const CGVector& globalVelocity, int cgi, int cgShift); +template +Eigen::Matrix cgLocal(const CGVector& globalVelocity, int cgi, + int cgShift); -template <> +template<> Eigen::Matrix cgLocal(const CGVector<1>& vGlobal, int cgi, int cgShift) { Eigen::Matrix vLocal; @@ -322,17 +317,17 @@ Eigen::Matrix cgLocal(const CGVector<1>& vGlobal, int cgi, return vLocal; } -template <> +template<> Eigen::Matrix cgLocal(const CGVector<2>& vGlobal, int cgi, int cgShift) { Eigen::Matrix vLocal; - vLocal << vGlobal(cgi), vGlobal(cgi + 1), vGlobal(cgi + 2), vGlobal(cgi + cgShift), - vGlobal(cgi + 1 + cgShift), vGlobal(cgi + 2 + cgShift), vGlobal(cgi + 2 * cgShift), - vGlobal(cgi + 1 + 2 * cgShift), vGlobal(cgi + 2 + 2 * cgShift); + vLocal << vGlobal(cgi), vGlobal(cgi + 1), vGlobal(cgi + 2), vGlobal(cgi + cgShift), vGlobal( + cgi + 1 + cgShift), vGlobal(cgi + 2 + cgShift), vGlobal(cgi + 2 * cgShift), vGlobal( + cgi + 1 + 2 * cgShift), vGlobal(cgi + 2 + 2 * cgShift); return vLocal; } -template void CGDynamicsKernel::projectVelocityToStrain() +template void CGDynamicsKernel::projectVelocityToStrain() { // !!! must still be converted to the spherical system!!! @@ -344,17 +339,16 @@ template void CGDynamicsKernel::projectVelocityTo int dgi = smesh->nx * row; //!< Index of dg vector int cgi = CGdegree * cgshift * row; //!< Lower left index of cg vector - for (size_t col = 0; col < smesh->nx; - ++col, ++dgi, cgi += CGdegree) { // loop over all elements + for (size_t col = 0; col < smesh->nx; ++col, ++dgi, cgi += CGdegree) { // loop over all elements if (smesh->landmask[dgi] == 0) // only on ice continue; // get the local x/y - velocity coefficients on the element - Eigen::Matrix vx_local - = cgLocal(u, cgi, cgshift); - Eigen::Matrix vy_local - = cgLocal(v, cgi, cgshift); + Eigen::Matrix vx_local = cgLocal(u, cgi, + cgshift); + Eigen::Matrix vy_local = cgLocal(v, cgi, + cgshift); // Solve (E, Psi) = (0.5(DV + DV^T), Psi) // by integrating rhs and inverting with dG(stress) mass matrix @@ -371,22 +365,21 @@ template void CGDynamicsKernel::projectVelocityTo } } -template -void CGDynamicsKernel::addStressTensorCell( - const size_t eid, const size_t cx, const size_t cy) +template +void CGDynamicsKernel::addStressTensorCell(const size_t eid, const size_t cx, + const size_t cy) { Eigen::Vector tx = (pmap->divS1[eid] * s11.row(eid).transpose() - + pmap->divS2[eid] * s12.row(eid).transpose()); + + pmap->divS2[eid] * s12.row(eid).transpose()); Eigen::Vector ty = (pmap->divS1[eid] * s12.row(eid).transpose() - + pmap->divS2[eid] * s22.row(eid).transpose()); + + pmap->divS2[eid] * s22.row(eid).transpose()); if (smesh->CoordinateSystem == SPHERICAL) { tx += pmap->divM[eid] * s12.row(eid).transpose(); ty -= pmap->divM[eid] * s11.row(eid).transpose(); } const size_t cgRow = CGdegree * smesh->nx + 1; - const size_t cg_i - = CGdegree * cgRow * cy + CGdegree * cx; //!< lower left CG-index in element (cx,cy) + const size_t cg_i = CGdegree * cgRow * cy + CGdegree * cx; //!< lower left CG-index in element (cx,cy) // Fill the stress divergence values for (size_t row = 0; row <= CGdegree; ++row) { @@ -397,11 +390,11 @@ void CGDynamicsKernel::addStressTensorCell( } } -template void CGDynamicsKernel::stressDivergence() +template void CGDynamicsKernel::stressDivergence() { // Somewhat meaningless, but it uses the name in the former version of the code - auto& tx = dStressX; - auto& ty = dStressY; + auto &tx = dStressX; + auto &ty = dStressY; // Zero the stress gradient vectors #pragma omp parallel for @@ -432,7 +425,7 @@ template void CGDynamicsKernel::stressDivergence( VectorManipulations::CGAveragePeriodic(*smesh, ty); } -template +template void CGDynamicsKernel::dirichletZero(CGVector& v) const { // TR 07.04.2025: Dirichlet Zero (u=v=0) holds on land and on the boundary betreen @@ -444,16 +437,16 @@ void CGDynamicsKernel::dirichletZero(CGVector& v) const v(i) = 0; } -template void CGDynamicsKernel::applyBoundaries() +template void CGDynamicsKernel::applyBoundaries() { dirichletZero(u); dirichletZero(v); // TODO Periodic boundary conditions. } -template -DGVector& CGDynamicsKernel::advectDGVField( - double timestep, DGVector& field, double lowerLimit, double upperLimit) +template +DGVector& CGDynamicsKernel::advectDGVField(double timestep, + DGVector& field, double lowerLimit, double upperLimit) { dgtransport->step(timestep, field); if (lowerLimit > -std::numeric_limits::infinity()) @@ -464,9 +457,9 @@ DGVector& CGDynamicsKernel::advectDGVField( return field; } -template -void CGDynamicsKernel::updateIceOceanStress( - const CGVector& uIce, const CGVector& vIce) +template +void CGDynamicsKernel::updateIceOceanStress(const CGVector& uIce, + const CGVector& vIce) { const double FOcean = baseParams.COcean * baseParams.rhoOcean; @@ -481,8 +474,8 @@ void CGDynamicsKernel::updateIceOceanStress( } // Instantiate the templates for all (1, 3, 6) degrees of DGadvection -template class CGDynamicsKernel<1>; -template class CGDynamicsKernel<3>; -template class CGDynamicsKernel<6>; +template class CGDynamicsKernel<1> ; +template class CGDynamicsKernel<3> ; +template class CGDynamicsKernel<6> ; } diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index 7aaca18d9..970ff3745 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -18,7 +18,7 @@ namespace Nextsim { // The brittle momentum solver for CG velocity fields -template class BrittleCGDynamicsKernel : public CGDynamicsKernel { +template class BrittleCGDynamicsKernel: public CGDynamicsKernel { protected: using DynamicsKernel::hice; using DynamicsKernel::cice; @@ -58,12 +58,10 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern public: using CGDynamicsKernel::advectField; - BrittleCGDynamicsKernel( - StressUpdateStep& stressStepIn, const BBMParameters& paramsIn) - : CGDynamicsKernel(paramsIn) - , stressStep(stressStepIn) - , params(paramsIn) - , stresstransport(nullptr) + BrittleCGDynamicsKernel(StressUpdateStep& stressStepIn, + const BBMParameters& paramsIn) + : CGDynamicsKernel(paramsIn), stressStep(stressStepIn), params(paramsIn), stresstransport( + nullptr) { } virtual ~BrittleCGDynamicsKernel() = default; @@ -85,7 +83,10 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern } // The brittle rheologies use avgU and avgV to do the advection, not u and v, like mEVP - void prepareAdvection() override { dgtransport->prepareAdvection(avgU, avgV); } + void prepareAdvection() override + { + dgtransport->prepareAdvection(avgU, avgV); + } /*! * Advection function for stress. Stress components have no limits. @@ -116,7 +117,7 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern { advectDynamicsFields(tst.step.seconds()); - prepareIteration({ { hiceName, hice }, { ciceName, cice } }); + prepareIteration( { { hiceName, hice }, { ciceName, cice } }); // The timestep for the brittle solver is the solver subtimestep deltaT = tst.step.seconds() / params.nSteps; @@ -128,11 +129,11 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern projectVelocityToStrain(); - std::array>, N_TENSOR_ELEMENTS> stress - = { s11, s12, s22 }; // Call the step function on the StressUpdateStep class + std::array>, N_TENSOR_ELEMENTS> stress = { + s11, s12, s22 }; // Call the step function on the StressUpdateStep class // Call the step function on the StressUpdateStep class - stressStep.stressUpdateHighOrder( - params, *smesh, stress, { e11, e12, e22 }, hice, cice, deltaT); + stressStep.stressUpdateHighOrder(params, *smesh, stress, { e11, e12, e22 }, hice, cice, + deltaT); stressDivergence(); // Compute divergence of stress tensor @@ -177,8 +178,8 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern CGVector avgU; CGVector avgV; - StressUpdateStep& stressStep; - const BBMParameters& params; + StressUpdateStep &stressStep; + const BBMParameters ¶ms; std::unique_ptr> stresstransport; @@ -216,22 +217,22 @@ template class BrittleCGDynamicsKernel : public CGDynamicsKern // Atmospheric drag const double dragAtm = cgA(i) * FAtm * std::hypot(uAtmos(i), vAtmos(i)); const double tauX = dragAtm * uAtmos(i) - + cPrime * (uOcean(i) * cosOceanAngle - vOcean(i) * sinOceanAngle); + + cPrime * (uOcean(i) * cosOceanAngle - vOcean(i) * sinOceanAngle); const double tauY = dragAtm * vAtmos(i) - + cPrime * (vOcean(i) * cosOceanAngle + uOcean(i) * sinOceanAngle); + + cPrime * (vOcean(i) * cosOceanAngle + uOcean(i) * sinOceanAngle); // Stress gradient const double gradX = dStressX(i) / pmap->lumpedcgmass(i) - - params.rhoIce * cgH(i) * PhysicalConstants::g * xGradSeaSurfaceHeight(i); + - params.rhoIce * cgH(i) * PhysicalConstants::g * xGradSeaSurfaceHeight(i); const double gradY = dStressY(i) / pmap->lumpedcgmass(i) - - params.rhoIce * cgH(i) * PhysicalConstants::g * yGradSeaSurfaceHeight(i); + - params.rhoIce * cgH(i) * PhysicalConstants::g * yGradSeaSurfaceHeight(i); u(i) = alpha * uIce + beta * vIce - + dteOverMass * (alpha * (gradX + tauX) + beta * (gradY + tauY)); + + dteOverMass * (alpha * (gradX + tauX) + beta * (gradY + tauY)); u(i) *= rDenom; v(i) = alpha * vIce - beta * uIce - + dteOverMass * (alpha * (gradY + tauY) + beta * (gradX + tauX)); + + dteOverMass * (alpha * (gradY + tauY) + beta * (gradX + tauX)); v(i) *= rDenom; // Calculate the contribution to the average velocity diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index 21d371c8f..9d51048ae 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -31,9 +31,9 @@ namespace Nextsim { // forward declare the class holding the potentially non-DG parts -template class DynamicsInternals; +template class DynamicsInternals; -template class DynamicsKernel { +template class DynamicsKernel { public: typedef std::pair&> DataMapping; typedef std::map DataMap; @@ -44,7 +44,7 @@ template class DynamicsKernel { { //! Define the spatial mesh smesh = std::make_unique( - (isSpherical) ? Nextsim::SPHERICAL : Nextsim::CARTESIAN); + (isSpherical) ? Nextsim::SPHERICAL : Nextsim::CARTESIAN); smesh->coordinatesFromModelArray(coords); if (isSpherical) @@ -64,7 +64,7 @@ template class DynamicsKernel { // s12.resize_by_mesh(*smesh); // s22.resize_by_mesh(*smesh); - // Set initial values to zero. Prognostic fields will be filled from the restart file. +// Set initial values to zero. Prognostic fields will be filled from the restart file. // e11.zero(); // e12.zero(); // e22.zero(); @@ -91,14 +91,8 @@ template class DynamicsKernel { */ virtual void setData(const std::string& name, const ModelArray& data) { - static const std::set vectorHolderNames = { - hiceName, - ciceName, - hsnowName, - stress11Name, - stress12Name, - stress22Name, - }; + static const std::set vectorHolderNames = { hiceName, ciceName, hsnowName, + stress11Name, stress12Name, stress22Name, }; // Special cases: hice, cice, (damage, stress) <- not yet implemented if (vectorHolderNames.count(name)) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); @@ -131,14 +125,17 @@ template class DynamicsKernel { HField data(ModelArray::Type::H); if (name == hiceName || name == ciceName || name == hsnowName) { throw std::runtime_error( - std::string("DynamicsKernel::getDG0Data: Use array sharing for ") + name); + std::string("DynamicsKernel::getDG0Data: Use array sharing for ") + name); } else { // Any other named field must exist return ModelArray(ModelArray::component0Type(ModelArray::AdvectionType)); } } - virtual void update(const TimestepTime& tst) { ++stepNumber; } + virtual void update(const TimestepTime& tst) + { + ++stepNumber; + } /*! * Prepares the transport objects to perform the advection step. @@ -156,9 +153,9 @@ template class DynamicsKernel { * * @return A reference to the advected array. */ - ModelArray& advectField(double timestep, ModelArray& field, - double lowerLimit = -std::numeric_limits::infinity(), - double upperLimit = std::numeric_limits::infinity()) + ModelArray& advectField(double timestep, ModelArray& field, double lowerLimit = + -std::numeric_limits::infinity(), double upperLimit = + std::numeric_limits::infinity()) { DGVectorHolder holder(field); advectDGVField(timestep, holder, lowerLimit, upperLimit); @@ -176,9 +173,9 @@ template class DynamicsKernel { * @return A reference to the advected array. */ virtual DGVector& advectDGVField(double timestep, DGVector& field, - double lowerLimit = -std::numeric_limits::infinity(), - double upperLimit = std::numeric_limits::infinity()) - = 0; + double lowerLimit = -std::numeric_limits::infinity(), double upperLimit = + std::numeric_limits::infinity()) + = 0; virtual void advectDynamicsFields(double timestep) { diff --git a/dynamics/src/include/dgVectorHolder.hpp b/dynamics/src/include/dgVectorHolder.hpp index b41585aeb..ac6c9f2ad 100644 --- a/dynamics/src/include/dgVectorHolder.hpp +++ b/dynamics/src/include/dgVectorHolder.hpp @@ -12,46 +12,66 @@ namespace Nextsim { -template class DGVectorHolder { +template class DGVectorHolder { public: using EigenDGVector = typename DGVector::EigenDGVector; DGVectorHolder() - : ref(nullptr) + : ref(nullptr) { } DGVectorHolder(ModelArray& ma) - : DGVectorHolder(static_cast(ma)) + : DGVectorHolder(static_cast(ma)) { if (ma.nComponents() != DG) { - std::string msg = "DGVectorHolder<"+std::to_string(DG)+">(ModelArray&): incorrect number of components = " + std::to_string(ma.nComponents()); + std::string msg = "DGVectorHolder<" + std::to_string(DG) + + ">(ModelArray&): incorrect number of components = " + + std::to_string(ma.nComponents()); throw std::length_error(msg); } } DGVectorHolder(ModelArray::DataType& madt) - : DGVectorHolder(reinterpret_cast(madt)) + : DGVectorHolder(reinterpret_cast(madt)) { } DGVectorHolder(EigenDGVector& edgv) - : ref(&edgv) + : ref(&edgv) { } DGVectorHolder(DGVector& dgv) - : ref(&dgv) + : ref(&dgv) { } - operator DGVector&() { return reinterpret_cast&>(*ref); } - operator const DGVector&() const { return reinterpret_cast&>(*ref); } + operator DGVector&() + { + return reinterpret_cast&>(*ref); + } + operator const DGVector&() const + { + return reinterpret_cast&>(*ref); + } - double& operator()(size_t i, size_t j) { return (*ref)(i, j); } + double& operator()(size_t i, size_t j) + { + return (*ref)(i, j); + } - void zero() { ref->setZero(); } + void zero() + { + ref->setZero(); + } - auto row(Eigen::Index i) { return static_cast&>(*ref).row(i); } - auto col(Eigen::Index i) { return static_cast&>(*ref).col(i); } + auto row(Eigen::Index i) + { + return static_cast&>(*ref).row(i); + } + auto col(Eigen::Index i) + { + return static_cast&>(*ref).col(i); + } private: - EigenDGVector* ref; + EigenDGVector *ref; }; } From 0707c10b8b9ac160445d2e26075d5f2da8e94c49 Mon Sep 17 00:00:00 2001 From: Tim Spain Date: Wed, 24 Sep 2025 09:23:11 +0200 Subject: [PATCH 5/5] clang formatting --- dynamics/src/CGDynamicsKernel.cpp | 148 +++++++++--------- .../src/include/BrittleCGDynamicsKernel.hpp | 55 ++++--- dynamics/src/include/DynamicsKernel.hpp | 69 ++++---- dynamics/src/include/dgVectorHolder.hpp | 48 ++---- 4 files changed, 156 insertions(+), 164 deletions(-) diff --git a/dynamics/src/CGDynamicsKernel.cpp b/dynamics/src/CGDynamicsKernel.cpp index b28dc58d8..94e4a46dc 100644 --- a/dynamics/src/CGDynamicsKernel.cpp +++ b/dynamics/src/CGDynamicsKernel.cpp @@ -11,25 +11,25 @@ #include "include/ModelArray.hpp" #include "include/constants.hpp" +#include "include/DGModelArray.hpp" #include "include/Interpolations.hpp" #include "include/ParametricMap.hpp" #include "include/VectorManipulations.hpp" #include "include/cgVector.hpp" -#include "include/DGModelArray.hpp" #include namespace Nextsim { -template +template CGDynamicsKernel::CGDynamicsKernel(const DynamicsParameters& params) - : baseParams(params) + : baseParams(params) { } -template -void CGDynamicsKernel::initialise(const ModelArray& coords, bool isSpherical, - const ModelArray& mask) +template +void CGDynamicsKernel::initialise( + const ModelArray& coords, bool isSpherical, const ModelArray& mask) { DynamicsKernel::initialise(coords, isSpherical, mask); @@ -68,12 +68,17 @@ void CGDynamicsKernel::initialise(const ModelArray& coords, bool is sinOceanAngle = std::sin(radians(baseParams.oceanTurningAngle)); } -template +template void CGDynamicsKernel::setData(const std::string& name, const ModelArray& data) { - const std::map*> targetMap = { { uName, &u }, { vName, &v }, { - uWindName, &uAtmos }, { vWindName, &vAtmos }, { uOceanName, &uOcean }, { vOceanName, - &vOcean }, }; + const std::map*> targetMap = { + { uName, &u }, + { vName, &v }, + { uWindName, &uAtmos }, + { vWindName, &vAtmos }, + { uOceanName, &uOcean }, + { vOceanName, &vOcean }, + }; if (targetMap.count(name)) { ma2cg(data, *targetMap.at(name)); } else { @@ -81,9 +86,9 @@ void CGDynamicsKernel::setData(const std::string& name, const Model } } -template -void CGDynamicsKernel::setDGArray(const std::string& name, - ModelArray::DataType& dgData) +template +void CGDynamicsKernel::setDGArray( + const std::string& name, ModelArray::DataType& dgData) { if (name == stress11Name) { s11 = DGVectorHolder(dgData); @@ -96,7 +101,7 @@ void CGDynamicsKernel::setDGArray(const std::string& name, } } -template +template ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) const { if (name == uName) { @@ -128,28 +133,30 @@ ModelArray CGDynamicsKernel::getDG0Data(const std::string& name) co } else if (name == sigmaIName) { ModelArray data(ModelArray::Type::H); return DGModelArray::dg2ma( - Tools::TensorInvI(*smesh, static_cast&>(s11), - static_cast&>(s12), - static_cast&>(s22)), data); + Tools::TensorInvI(*smesh, static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22)), + data); } else if (name == sigmaIIName) { ModelArray data(ModelArray::Type::H); return DGModelArray::dg2ma( - Tools::TensorInvII(*smesh, static_cast&>(s11), - static_cast&>(s12), - static_cast&>(s22)), data); + Tools::TensorInvII(*smesh, static_cast&>(s11), + static_cast&>(s12), + static_cast&>(s22)), + data); } else { return DynamicsKernel::getDG0Data(name); } } -template void CGDynamicsKernel::prepareAdvection() +template void CGDynamicsKernel::prepareAdvection() { dgtransport->prepareAdvection(u, v); } -template +template void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( - const DGVector<1>& seaSurfaceHeight) + const DGVector<1>& seaSurfaceHeight) { // First transform to CG1 Vector and do all computations in CG1 CGVector<1> cgSeasurfaceHeight(*smesh); @@ -171,9 +178,8 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( for (size_t cx = 0; cx < smesh->nx; ++cx, ++eid, ++cg1id) { // get local CG nodes Eigen::Vector loc_cgSSH = { cgSeasurfaceHeight(cg1id), - cgSeasurfaceHeight(cg1id + 1), cgSeasurfaceHeight( - cg1id + smesh->nx + 1), cgSeasurfaceHeight( - cg1id + smesh->nx + 1 + 1) }; + cgSeasurfaceHeight(cg1id + 1), cgSeasurfaceHeight(cg1id + smesh->nx + 1), + cgSeasurfaceHeight(cg1id + smesh->nx + 1 + 1) }; // compute grad Eigen::Vector tx = pmap->dX_SSH[eid] * loc_cgSSH; @@ -204,14 +210,14 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( size_t cg1row = smesh->nx + 1; size_t topleft = smesh->ny * cg1row; for (size_t i = 1; i < smesh->nx; ++i) // bottom / top - { + { uGrad(i) = uGrad(i + cg1row); vGrad(i) = vGrad(i + cg1row); uGrad(topleft + i) = uGrad(topleft + i - cg1row); vGrad(topleft + i) = vGrad(topleft + i - cg1row); } for (size_t i = 1; i < smesh->ny; ++i) // left / right - { + { uGrad(i * cg1row) = uGrad(i * cg1row + 1); vGrad(i * cg1row) = vGrad(i * cg1row + 1); @@ -246,7 +252,7 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( // along lines #pragma omp parallel for for (size_t iy = 0; iy <= smesh->ny; ++iy) // horizontal - { + { size_t icg1 = (smesh->nx + 1) * iy; size_t icg2 = (2 * smesh->nx + 1) * 2 * iy + 1; for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) { @@ -256,7 +262,7 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( } #pragma omp parallel for for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical - { + { size_t icg1 = (smesh->nx + 1) * iy; size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1); for (size_t ix = 0; ix <= smesh->nx; ++ix, ++icg1, icg2 += 2) { @@ -268,22 +274,22 @@ void CGDynamicsKernel::computeGradientOfSeaSurfaceHeight( // midpoints #pragma omp parallel for for (size_t iy = 0; iy < smesh->ny; ++iy) // vertical - { + { size_t icg1 = (smesh->nx + 1) * iy; size_t icg2 = (2 * smesh->nx + 1) * (2 * iy + 1) + 1; for (size_t ix = 0; ix < smesh->nx; ++ix, ++icg1, icg2 += 2) { xGradSeaSurfaceHeight(icg2) = 0.25 - * (uGrad(icg1) + uGrad(icg1 + 1) + uGrad(icg1 + cg1row) - + uGrad(icg1 + cg1row + 1)); + * (uGrad(icg1) + uGrad(icg1 + 1) + uGrad(icg1 + cg1row) + + uGrad(icg1 + cg1row + 1)); yGradSeaSurfaceHeight(icg2) = 0.25 - * (vGrad(icg1) + vGrad(icg1 + 1) + vGrad(icg1 + cg1row) - + vGrad(icg1 + cg1row + 1)); + * (vGrad(icg1) + vGrad(icg1 + 1) + vGrad(icg1 + cg1row) + + vGrad(icg1 + cg1row + 1)); } } } } -template void CGDynamicsKernel::prepareIteration(const DataMap& data) +template void CGDynamicsKernel::prepareIteration(const DataMap& data) { // interpolate ice height and concentration to local cg Variables Interpolations::DG2CG(*smesh, cgH, data.at(hiceName)); @@ -305,11 +311,11 @@ template void CGDynamicsKernel::prepareIteration(c cgH = cgH.cwiseMax(0.05); } -template -Eigen::Matrix cgLocal(const CGVector& globalVelocity, int cgi, - int cgShift); +template +Eigen::Matrix cgLocal( + const CGVector& globalVelocity, int cgi, int cgShift); -template<> +template <> Eigen::Matrix cgLocal(const CGVector<1>& vGlobal, int cgi, int cgShift) { Eigen::Matrix vLocal; @@ -317,17 +323,17 @@ Eigen::Matrix cgLocal(const CGVector<1>& vGlobal, int cgi, return vLocal; } -template<> +template <> Eigen::Matrix cgLocal(const CGVector<2>& vGlobal, int cgi, int cgShift) { Eigen::Matrix vLocal; - vLocal << vGlobal(cgi), vGlobal(cgi + 1), vGlobal(cgi + 2), vGlobal(cgi + cgShift), vGlobal( - cgi + 1 + cgShift), vGlobal(cgi + 2 + cgShift), vGlobal(cgi + 2 * cgShift), vGlobal( - cgi + 1 + 2 * cgShift), vGlobal(cgi + 2 + 2 * cgShift); + vLocal << vGlobal(cgi), vGlobal(cgi + 1), vGlobal(cgi + 2), vGlobal(cgi + cgShift), + vGlobal(cgi + 1 + cgShift), vGlobal(cgi + 2 + cgShift), vGlobal(cgi + 2 * cgShift), + vGlobal(cgi + 1 + 2 * cgShift), vGlobal(cgi + 2 + 2 * cgShift); return vLocal; } -template void CGDynamicsKernel::projectVelocityToStrain() +template void CGDynamicsKernel::projectVelocityToStrain() { // !!! must still be converted to the spherical system!!! @@ -339,16 +345,17 @@ template void CGDynamicsKernel::projectVelocityToS int dgi = smesh->nx * row; //!< Index of dg vector int cgi = CGdegree * cgshift * row; //!< Lower left index of cg vector - for (size_t col = 0; col < smesh->nx; ++col, ++dgi, cgi += CGdegree) { // loop over all elements + for (size_t col = 0; col < smesh->nx; + ++col, ++dgi, cgi += CGdegree) { // loop over all elements if (smesh->landmask[dgi] == 0) // only on ice continue; // get the local x/y - velocity coefficients on the element - Eigen::Matrix vx_local = cgLocal(u, cgi, - cgshift); - Eigen::Matrix vy_local = cgLocal(v, cgi, - cgshift); + Eigen::Matrix vx_local + = cgLocal(u, cgi, cgshift); + Eigen::Matrix vy_local + = cgLocal(v, cgi, cgshift); // Solve (E, Psi) = (0.5(DV + DV^T), Psi) // by integrating rhs and inverting with dG(stress) mass matrix @@ -365,21 +372,22 @@ template void CGDynamicsKernel::projectVelocityToS } } -template -void CGDynamicsKernel::addStressTensorCell(const size_t eid, const size_t cx, - const size_t cy) +template +void CGDynamicsKernel::addStressTensorCell( + const size_t eid, const size_t cx, const size_t cy) { Eigen::Vector tx = (pmap->divS1[eid] * s11.row(eid).transpose() - + pmap->divS2[eid] * s12.row(eid).transpose()); + + pmap->divS2[eid] * s12.row(eid).transpose()); Eigen::Vector ty = (pmap->divS1[eid] * s12.row(eid).transpose() - + pmap->divS2[eid] * s22.row(eid).transpose()); + + pmap->divS2[eid] * s22.row(eid).transpose()); if (smesh->CoordinateSystem == SPHERICAL) { tx += pmap->divM[eid] * s12.row(eid).transpose(); ty -= pmap->divM[eid] * s11.row(eid).transpose(); } const size_t cgRow = CGdegree * smesh->nx + 1; - const size_t cg_i = CGdegree * cgRow * cy + CGdegree * cx; //!< lower left CG-index in element (cx,cy) + const size_t cg_i + = CGdegree * cgRow * cy + CGdegree * cx; //!< lower left CG-index in element (cx,cy) // Fill the stress divergence values for (size_t row = 0; row <= CGdegree; ++row) { @@ -390,11 +398,11 @@ void CGDynamicsKernel::addStressTensorCell(const size_t eid, const } } -template void CGDynamicsKernel::stressDivergence() +template void CGDynamicsKernel::stressDivergence() { // Somewhat meaningless, but it uses the name in the former version of the code - auto &tx = dStressX; - auto &ty = dStressY; + auto& tx = dStressX; + auto& ty = dStressY; // Zero the stress gradient vectors #pragma omp parallel for @@ -425,7 +433,7 @@ template void CGDynamicsKernel::stressDivergence() VectorManipulations::CGAveragePeriodic(*smesh, ty); } -template +template void CGDynamicsKernel::dirichletZero(CGVector& v) const { // TR 07.04.2025: Dirichlet Zero (u=v=0) holds on land and on the boundary betreen @@ -437,16 +445,16 @@ void CGDynamicsKernel::dirichletZero(CGVector& v) const v(i) = 0; } -template void CGDynamicsKernel::applyBoundaries() +template void CGDynamicsKernel::applyBoundaries() { dirichletZero(u); dirichletZero(v); // TODO Periodic boundary conditions. } -template -DGVector& CGDynamicsKernel::advectDGVField(double timestep, - DGVector& field, double lowerLimit, double upperLimit) +template +DGVector& CGDynamicsKernel::advectDGVField( + double timestep, DGVector& field, double lowerLimit, double upperLimit) { dgtransport->step(timestep, field); if (lowerLimit > -std::numeric_limits::infinity()) @@ -457,9 +465,9 @@ DGVector& CGDynamicsKernel::advectDGVField(double time return field; } -template -void CGDynamicsKernel::updateIceOceanStress(const CGVector& uIce, - const CGVector& vIce) +template +void CGDynamicsKernel::updateIceOceanStress( + const CGVector& uIce, const CGVector& vIce) { const double FOcean = baseParams.COcean * baseParams.rhoOcean; @@ -474,8 +482,8 @@ void CGDynamicsKernel::updateIceOceanStress(const CGVector ; -template class CGDynamicsKernel<3> ; -template class CGDynamicsKernel<6> ; +template class CGDynamicsKernel<1>; +template class CGDynamicsKernel<3>; +template class CGDynamicsKernel<6>; } diff --git a/dynamics/src/include/BrittleCGDynamicsKernel.hpp b/dynamics/src/include/BrittleCGDynamicsKernel.hpp index 970ff3745..eda40f28f 100644 --- a/dynamics/src/include/BrittleCGDynamicsKernel.hpp +++ b/dynamics/src/include/BrittleCGDynamicsKernel.hpp @@ -18,7 +18,7 @@ namespace Nextsim { // The brittle momentum solver for CG velocity fields -template class BrittleCGDynamicsKernel: public CGDynamicsKernel { +template class BrittleCGDynamicsKernel : public CGDynamicsKernel { protected: using DynamicsKernel::hice; using DynamicsKernel::cice; @@ -58,10 +58,12 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel public: using CGDynamicsKernel::advectField; - BrittleCGDynamicsKernel(StressUpdateStep& stressStepIn, - const BBMParameters& paramsIn) - : CGDynamicsKernel(paramsIn), stressStep(stressStepIn), params(paramsIn), stresstransport( - nullptr) + BrittleCGDynamicsKernel( + StressUpdateStep& stressStepIn, const BBMParameters& paramsIn) + : CGDynamicsKernel(paramsIn) + , stressStep(stressStepIn) + , params(paramsIn) + , stresstransport(nullptr) { } virtual ~BrittleCGDynamicsKernel() = default; @@ -83,10 +85,7 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel } // The brittle rheologies use avgU and avgV to do the advection, not u and v, like mEVP - void prepareAdvection() override - { - dgtransport->prepareAdvection(avgU, avgV); - } + void prepareAdvection() override { dgtransport->prepareAdvection(avgU, avgV); } /*! * Advection function for stress. Stress components have no limits. @@ -117,7 +116,7 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel { advectDynamicsFields(tst.step.seconds()); - prepareIteration( { { hiceName, hice }, { ciceName, cice } }); + prepareIteration({ { hiceName, hice }, { ciceName, cice } }); // The timestep for the brittle solver is the solver subtimestep deltaT = tst.step.seconds() / params.nSteps; @@ -129,11 +128,11 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel projectVelocityToStrain(); - std::array>, N_TENSOR_ELEMENTS> stress = { - s11, s12, s22 }; // Call the step function on the StressUpdateStep class + std::array>, N_TENSOR_ELEMENTS> stress + = { s11, s12, s22 }; // Call the step function on the StressUpdateStep class // Call the step function on the StressUpdateStep class - stressStep.stressUpdateHighOrder(params, *smesh, stress, { e11, e12, e22 }, hice, cice, - deltaT); + stressStep.stressUpdateHighOrder( + params, *smesh, stress, { e11, e12, e22 }, hice, cice, deltaT); stressDivergence(); // Compute divergence of stress tensor @@ -154,12 +153,12 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel { if (name == damageName) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); -// } else if (name == stress11Name) { -// DGModelArray::ma2dg(data, s11); -// } else if (name == stress12Name) { -// DGModelArray::ma2dg(data, s12); -// } else if (name == stress22Name) { -// DGModelArray::ma2dg(data, s22); + // } else if (name == stress11Name) { + // DGModelArray::ma2dg(data, s11); + // } else if (name == stress12Name) { + // DGModelArray::ma2dg(data, s12); + // } else if (name == stress22Name) { + // DGModelArray::ma2dg(data, s22); } else { CGDynamicsKernel::setData(name, data); } @@ -178,8 +177,8 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel CGVector avgU; CGVector avgV; - StressUpdateStep &stressStep; - const BBMParameters ¶ms; + StressUpdateStep& stressStep; + const BBMParameters& params; std::unique_ptr> stresstransport; @@ -217,22 +216,22 @@ template class BrittleCGDynamicsKernel: public CGDynamicsKernel // Atmospheric drag const double dragAtm = cgA(i) * FAtm * std::hypot(uAtmos(i), vAtmos(i)); const double tauX = dragAtm * uAtmos(i) - + cPrime * (uOcean(i) * cosOceanAngle - vOcean(i) * sinOceanAngle); + + cPrime * (uOcean(i) * cosOceanAngle - vOcean(i) * sinOceanAngle); const double tauY = dragAtm * vAtmos(i) - + cPrime * (vOcean(i) * cosOceanAngle + uOcean(i) * sinOceanAngle); + + cPrime * (vOcean(i) * cosOceanAngle + uOcean(i) * sinOceanAngle); // Stress gradient const double gradX = dStressX(i) / pmap->lumpedcgmass(i) - - params.rhoIce * cgH(i) * PhysicalConstants::g * xGradSeaSurfaceHeight(i); + - params.rhoIce * cgH(i) * PhysicalConstants::g * xGradSeaSurfaceHeight(i); const double gradY = dStressY(i) / pmap->lumpedcgmass(i) - - params.rhoIce * cgH(i) * PhysicalConstants::g * yGradSeaSurfaceHeight(i); + - params.rhoIce * cgH(i) * PhysicalConstants::g * yGradSeaSurfaceHeight(i); u(i) = alpha * uIce + beta * vIce - + dteOverMass * (alpha * (gradX + tauX) + beta * (gradY + tauY)); + + dteOverMass * (alpha * (gradX + tauX) + beta * (gradY + tauY)); u(i) *= rDenom; v(i) = alpha * vIce - beta * uIce - + dteOverMass * (alpha * (gradY + tauY) + beta * (gradX + tauX)); + + dteOverMass * (alpha * (gradY + tauY) + beta * (gradX + tauX)); v(i) *= rDenom; // Calculate the contribution to the average velocity diff --git a/dynamics/src/include/DynamicsKernel.hpp b/dynamics/src/include/DynamicsKernel.hpp index 9d51048ae..8271b4aad 100644 --- a/dynamics/src/include/DynamicsKernel.hpp +++ b/dynamics/src/include/DynamicsKernel.hpp @@ -31,9 +31,9 @@ namespace Nextsim { // forward declare the class holding the potentially non-DG parts -template class DynamicsInternals; +template class DynamicsInternals; -template class DynamicsKernel { +template class DynamicsKernel { public: typedef std::pair&> DataMapping; typedef std::map DataMap; @@ -44,7 +44,7 @@ template class DynamicsKernel { { //! Define the spatial mesh smesh = std::make_unique( - (isSpherical) ? Nextsim::SPHERICAL : Nextsim::CARTESIAN); + (isSpherical) ? Nextsim::SPHERICAL : Nextsim::CARTESIAN); smesh->coordinatesFromModelArray(coords); if (isSpherical) @@ -57,20 +57,20 @@ template class DynamicsKernel { seaSurfaceHeight.resize_by_mesh(*smesh); -// e11.resize_by_mesh(*smesh); -// e12.resize_by_mesh(*smesh); -// e22.resize_by_mesh(*smesh); -// s11.resize_by_mesh(*smesh); -// s12.resize_by_mesh(*smesh); -// s22.resize_by_mesh(*smesh); - -// Set initial values to zero. Prognostic fields will be filled from the restart file. -// e11.zero(); -// e12.zero(); -// e22.zero(); -// s11.zero(); -// s12.zero(); -// s22.zero(); + // e11.resize_by_mesh(*smesh); + // e12.resize_by_mesh(*smesh); + // e22.resize_by_mesh(*smesh); + // s11.resize_by_mesh(*smesh); + // s12.resize_by_mesh(*smesh); + // s22.resize_by_mesh(*smesh); + + // Set initial values to zero. Prognostic fields will be filled from the restart file. + // e11.zero(); + // e12.zero(); + // e22.zero(); + // s11.zero(); + // s12.zero(); + // s22.zero(); } /*! @@ -91,8 +91,14 @@ template class DynamicsKernel { */ virtual void setData(const std::string& name, const ModelArray& data) { - static const std::set vectorHolderNames = { hiceName, ciceName, hsnowName, - stress11Name, stress12Name, stress22Name, }; + static const std::set vectorHolderNames = { + hiceName, + ciceName, + hsnowName, + stress11Name, + stress12Name, + stress22Name, + }; // Special cases: hice, cice, (damage, stress) <- not yet implemented if (vectorHolderNames.count(name)) { throw std::runtime_error(std::string("Use setDGArray() to set the data for ") + name); @@ -125,17 +131,14 @@ template class DynamicsKernel { HField data(ModelArray::Type::H); if (name == hiceName || name == ciceName || name == hsnowName) { throw std::runtime_error( - std::string("DynamicsKernel::getDG0Data: Use array sharing for ") + name); + std::string("DynamicsKernel::getDG0Data: Use array sharing for ") + name); } else { // Any other named field must exist return ModelArray(ModelArray::component0Type(ModelArray::AdvectionType)); } } - virtual void update(const TimestepTime& tst) - { - ++stepNumber; - } + virtual void update(const TimestepTime& tst) { ++stepNumber; } /*! * Prepares the transport objects to perform the advection step. @@ -153,9 +156,9 @@ template class DynamicsKernel { * * @return A reference to the advected array. */ - ModelArray& advectField(double timestep, ModelArray& field, double lowerLimit = - -std::numeric_limits::infinity(), double upperLimit = - std::numeric_limits::infinity()) + ModelArray& advectField(double timestep, ModelArray& field, + double lowerLimit = -std::numeric_limits::infinity(), + double upperLimit = std::numeric_limits::infinity()) { DGVectorHolder holder(field); advectDGVField(timestep, holder, lowerLimit, upperLimit); @@ -173,9 +176,9 @@ template class DynamicsKernel { * @return A reference to the advected array. */ virtual DGVector& advectDGVField(double timestep, DGVector& field, - double lowerLimit = -std::numeric_limits::infinity(), double upperLimit = - std::numeric_limits::infinity()) - = 0; + double lowerLimit = -std::numeric_limits::infinity(), + double upperLimit = std::numeric_limits::infinity()) + = 0; virtual void advectDynamicsFields(double timestep) { @@ -195,9 +198,9 @@ template class DynamicsKernel { DGVector<1> seaSurfaceHeight; //! Vectors storing strain and stress components -// DGVector e11, e12, e22; -// DGVector /*s11,*/ s12, s22; -// DGVectorHolder s11, s12, s22; + // DGVector e11, e12, e22; + // DGVector /*s11,*/ s12, s22; + // DGVectorHolder s11, s12, s22; size_t stepNumber = 0; diff --git a/dynamics/src/include/dgVectorHolder.hpp b/dynamics/src/include/dgVectorHolder.hpp index ac6c9f2ad..7b8285ac6 100644 --- a/dynamics/src/include/dgVectorHolder.hpp +++ b/dynamics/src/include/dgVectorHolder.hpp @@ -12,66 +12,48 @@ namespace Nextsim { -template class DGVectorHolder { +template class DGVectorHolder { public: using EigenDGVector = typename DGVector::EigenDGVector; DGVectorHolder() - : ref(nullptr) + : ref(nullptr) { } DGVectorHolder(ModelArray& ma) - : DGVectorHolder(static_cast(ma)) + : DGVectorHolder(static_cast(ma)) { if (ma.nComponents() != DG) { std::string msg = "DGVectorHolder<" + std::to_string(DG) - + ">(ModelArray&): incorrect number of components = " - + std::to_string(ma.nComponents()); + + ">(ModelArray&): incorrect number of components = " + + std::to_string(ma.nComponents()); throw std::length_error(msg); } } DGVectorHolder(ModelArray::DataType& madt) - : DGVectorHolder(reinterpret_cast(madt)) + : DGVectorHolder(reinterpret_cast(madt)) { } DGVectorHolder(EigenDGVector& edgv) - : ref(&edgv) + : ref(&edgv) { } DGVectorHolder(DGVector& dgv) - : ref(&dgv) + : ref(&dgv) { } - operator DGVector&() - { - return reinterpret_cast&>(*ref); - } - operator const DGVector&() const - { - return reinterpret_cast&>(*ref); - } + operator DGVector&() { return reinterpret_cast&>(*ref); } + operator const DGVector&() const { return reinterpret_cast&>(*ref); } - double& operator()(size_t i, size_t j) - { - return (*ref)(i, j); - } + double& operator()(size_t i, size_t j) { return (*ref)(i, j); } - void zero() - { - ref->setZero(); - } + void zero() { ref->setZero(); } - auto row(Eigen::Index i) - { - return static_cast&>(*ref).row(i); - } - auto col(Eigen::Index i) - { - return static_cast&>(*ref).col(i); - } + auto row(Eigen::Index i) { return static_cast&>(*ref).row(i); } + auto col(Eigen::Index i) { return static_cast&>(*ref).col(i); } private: - EigenDGVector *ref; + EigenDGVector* ref; }; }