[PL] GWF; Anisotropic hydraulic conductivity.

ufz · Mar 9, 2019 · 6bf8029 · 6bf8029
1 parent c578d6c
commit 6bf8029
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 2 deletions.
diff --git a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
@@ -53,7 +53,7 @@ std::unique_ptr<Process> createGroundwaterFlowProcess(
     auto& hydraulic_conductivity = findParameter<double>(
         config,
         //! \ogs_file_param_special{prj__processes__process__GROUNDWATER_FLOW__hydraulic_conductivity}
-        "hydraulic_conductivity", parameters, 1);
+        "hydraulic_conductivity", parameters, 0 /*arbitrary many components*/);
 
     DBUG("Use '%s' as hydraulic conductivity parameter.",
          hydraulic_conductivity.name.c_str());

diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -28,6 +28,34 @@ namespace GroundwaterFlow
 {
 const unsigned NUM_NODAL_DOF = 1;
 
+template <int GlobalDim>
+Eigen::Matrix<double, GlobalDim, GlobalDim> hydraulicConductivity(
+    std::vector<double> const& values)
+{
+    auto const size{values.size()};
+    if (size == 1)  // This is a duplicate but preferred case for GlobalDim==1.
+    {
+        return Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity() *
+               values[0];
+    }
+    if (size == GlobalDim)
+    {
+        return Eigen::Map<Eigen::Matrix<double, GlobalDim, 1> const>(
+                   values.data(), GlobalDim, 1)
+            .asDiagonal();
+    }
+    if (size == GlobalDim * GlobalDim)
+    {
+        return Eigen::Map<Eigen::Matrix<double, GlobalDim, GlobalDim> const>(
+            values.data(), GlobalDim, GlobalDim);
+    }
+
+    OGS_FATAL(
+        "Hydraulic conductivity parameter values size is neither one nor %d "
+        "nor %d squared, but %d.",
+        GlobalDim, GlobalDim, values.size());
+}
+
 class GroundwaterFlowLocalAssemblerInterface
     : public ProcessLib::LocalAssemblerInterface,
       public NumLib::ExtrapolatableElement
@@ -95,7 +123,8 @@ class LocalAssemblerData : public GroundwaterFlowLocalAssemblerInterface
             pos.setIntegrationPoint(ip);
             auto const& sm = _shape_matrices[ip];
             auto const& wp = _integration_method.getWeightedPoint(ip);
-            auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
+            auto const k = hydraulicConductivity<GlobalDim>(
+                _process_data.hydraulic_conductivity(t, pos));
 
             local_K.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ *
                                  sm.integralMeasure * wp.getWeight();