Some code refactoring to better manage Tetra (#72)

Author: danielemoretto44

Code/boundaries/Boundaries.m

@@ -8,6 +8,9 @@
 
   properties (Access = private)
     grid
+    bcList 
+    % to ensure neumann is applied before dirichlet in case of overlapping
+    % entity definition
   end
 
   methods (Access = public)
@@ -254,6 +257,13 @@ function removeBCentities(obj,bcId,list)
       obj.db(bcId) = bc;
 
     end
+
+
+    function bcList = getBCList(obj)
+
+      bcList = obj.bcList;
+
+    end
   end
 
 
@@ -276,7 +286,9 @@ function readInputFile(obj,inputStruct)
         inputStruct = readstruct(inputStruct.fileName,AttributeSuffix="");
       end
 
-      for i = 1:numel(inputStruct.BC)
+      nBC = numel(inputStruct.BC);
+
+      for i = 1:nBC
         % process each BC
         in = inputStruct.BC(i);
 
@@ -334,6 +346,16 @@ function readInputFile(obj,inputStruct)
         obj.db(name) = bc;
       end
 
+      % set correct order of boundary conditions. Dirichlet last
+      bcTypeList = [];
+      for bcId = string(obj.db.keys)
+        bcTypeList = [bcTypeList, obj.getType(bcId)];
+      end
+      idxDir  = strcmp(bcTypeList, 'Dirichlet'); 
+      bcOrd = [find(~idxDir), find(idxDir)];
+      bcNames = obj.db.keys;
+      obj.bcList = string(bcNames(bcOrd));
+
     end
 
   end

Code/boundaries/linkBoundSurf2TPFAFace.m

@@ -1,7 +1,7 @@
 function linkBoundSurf2TPFAFace(solver)
 % TO DO: call this from within the SinglePhaseFlow solver!
 
-bcs = solver.bcs;
+bcs = solver.domain.bcs;
 
 keys = bcs.db.keys;
 flRenum = false(length(keys),1);

Code/discretizer/Discretizer.m

@@ -1,7 +1,7 @@
 classdef Discretizer < handle
   % General discretizer class
 
-  properties (GetAccess=public, SetAccess=private)
+  properties (GetAccess=public, SetAccess=public)
     physicsSolvers                     % physics solvers database
     solverNames
     dofm
@@ -52,9 +52,9 @@
     end
 
     function applyBC(obj,t)
-      bcList = obj.bcs.db.keys;
+      bcList = obj.bcs.getBCList();
 
-      for bcId = string(bcList)
+      for bcId = bcList
         % loop over available bcs
         bcVar = obj.bcs.getVariable(bcId);
 
@@ -74,9 +74,9 @@ function applyBC(obj,t)
 
 
     function applyDirVal(obj,t)
-      bcList = obj.bcs.db.keys;
+      bcList = obj.bcs.getBCList;
 
-      for bcId = string(bcList)
+      for bcId = bcList
         % discard non-Dirichlet BC
         if ~strcmp(obj.bcs.getType(bcId),"Dirichlet")
           continue
@@ -458,12 +458,20 @@ function setInput(obj, varargin)
           case 'grid'
             obj.grid = value;
             % check that grid has been defined correctly
-            isGridCorrect = all([isfield(obj.grid,"topology");...
-              isfield(obj.grid,"cells");...
-              isfield(obj.grid,"faces")]);
+            if ~isfield(obj.grid,"topology")
+              obj.grid.topology = [];
+            end
+            if ~isfield(obj.grid,"cells")
+              obj.grid.cells = [];
+            end
+            if ~isfield(obj.grid,"faces")
+              obj.grid.faces = [];
+            end
+
+            isGridCorrect = numel(fieldnames(obj.grid))==3;
 
             assert(isGridCorrect,"Error in Discretizer: grid input is not correct. " + ...
-              "See the default value of the grid property in Discretizer.");
+              "grid must be a struct with fields: 'topology','cells','faces'.");
 
           case 'materials'
             assert(isa(value, 'Materials'),msg)

Code/discretizer/PhysicsSolver.m

@@ -33,27 +33,21 @@
   properties (GetAccess=public, SetAccess=protected)
     % handle to domain properties
     domain
-    dofm
+    simparams
     mesh
     elements
     faces
-    bcs
-    materials
-    simparams
   end
 
   methods
     function obj = PhysicsSolver(domain)
 
       % inputStruct: struct with additional solver-specific parameters
       obj.domain = domain;
-      obj.dofm = domain.dofm;
+      obj.simparams = domain.simparams;
       obj.mesh = domain.grid.topology;
-      obj.elements = domain.grid.cells;
       obj.faces = domain.grid.faces;
-      obj.materials = domain.materials;
-      obj.bcs = domain.bcs;
-      obj.simparams = domain.simparams;
+      obj.elements = domain.grid.cells;
 
     end
   end
@@ -165,8 +159,8 @@ function applyNeuBC(obj,bcId,bcDofs,bcVals)
 
       % Base application of Neumann boundary condition to the rhs.
       % bc values are subtracted since we solve du = J\(-rhs)
-      bcVar = obj.bcs.getVariable(bcId);
-      bcId = obj.dofm.getVariableId(bcVar);
+      bcVar = obj.domain.bcs.getVariable(bcId);
+      bcId = obj.domain.dofm.getVariableId(bcVar);
 
       % remove inactive dofs
       id = bcDofs == 0;
@@ -197,10 +191,10 @@ function applyDirBC(obj,bcId,bcDofs,bcVals)
       % sort bcDofs to improve sparse access performance
       [bcDofs,sortId] = sort(bcDofs);
 
-      bcVar = obj.bcs.getVariable(bcId);
-      bcVarId = obj.dofm.getVariableId(bcVar);
+      bcVar = obj.domain.bcs.getVariable(bcId);
+      bcVarId = obj.domain.dofm.getVariableId(bcVar);
 
-      nV = getNumberOfVariables(obj.dofm);
+      nV = getNumberOfVariables(obj.domain.dofm);
 
       % zero out rows (use transpose trick)
       for j = 1:nV
@@ -265,7 +259,7 @@ function applyDirBC(obj,bcId,bcDofs,bcVals)
 
     function out = BCapplies(obj,bcId)
 
-      bcVar = obj.bcs.getVariable(bcId);
+      bcVar = obj.domain.bcs.getVariable(bcId);
       out = any(strcmp(obj.getField(),bcVar));
 
     end

Code/discretizer/PhysicsSolvers/Coupling/BiotFullySaturated.m

@@ -30,10 +30,12 @@
 
     function registerSolver(obj,input)
 
+      dofm = obj.domain.dofm;
+
       % Register mechanics
       obj.mechSolver = Poromechanics(obj.domain);
       registerSolver(obj.mechSolver,input.(class(obj.mechSolver)));
-      obj.fldMech = obj.dofm.getVariableId(obj.mechSolver.getField());
+      obj.fldMech = dofm.getVariableId(obj.mechSolver.getField());
 
       % setup the solver with custom input
       if isfield(input,"SinglePhaseFlowFEM")
@@ -46,7 +48,7 @@ function registerSolver(obj,input)
       % Register fluids
       obj.flowScheme = obj.flowSolver.typeDiscretization();
       registerSolver(obj.flowSolver,input.(class(obj.flowSolver)));
-      obj.fldFlow = obj.dofm.getVariableId(obj.flowSolver.getField());
+      obj.fldFlow = dofm.getVariableId(obj.flowSolver.getField());
 
     end
 
@@ -85,8 +87,10 @@ function computeMatBiot(obj)
       % compute coupling matrix only where mechanics and flow are
       % active
 
-      cellTagFlow = obj.dofm.getTargetRegions(obj.fldFlow);
-      cellTagMech = obj.dofm.getTargetRegions(obj.fldMech);
+      dofm = obj.domain.dofm;
+
+      cellTagFlow = dofm.getTargetRegions(obj.fldFlow);
+      cellTagMech = dofm.getTargetRegions(obj.fldMech);
 
       % find cell tag where both flow and mechanics are active
       cellTags = intersect(cellTagMech,cellTagFlow);
@@ -103,14 +107,14 @@ function computeMatBiot(obj)
       end
 
       [iiVec,jjVec,Qvec] = deal(zeros(nEntries,1));
-      nDoF1 = obj.dofm.getNumbDoF(obj.fldMech);
-      nDoF2 = obj.dofm.getNumbDoF(obj.fldFlow);
+      nDoF1 = dofm.getNumbDoF(obj.fldMech);
+      nDoF2 = dofm.getNumbDoF(obj.fldFlow);
       % consider replacing the string field with an integer
 
       l1 = 0;
       for el=subCells'
-
-        biot = obj.materials.getMaterial(obj.mesh.cellTag(el)).PorousRock.getBiotCoefficient();
+        mat = obj.domain.materials.getMaterial(obj.mesh.cellTag(el));
+        biot = mat.PorousRock.getBiotCoefficient();
         elem = getElement(obj.elements,obj.mesh.cellVTKType(el));
         nG = elem.GaussPts.nNode;
         nodes = obj.mesh.cells(el,1:obj.mesh.cellNumVerts(el));
@@ -121,7 +125,7 @@ function computeMatBiot(obj)
         B = zeros(6,elem.nNode*obj.mesh.nDim,nG);
         B(elem.indB(:,2)) = N(elem.indB(:,1));
         Nref = getBasisFinGPoints(elem);
-        dofrow = getLocalDoF(obj.dofm,obj.fldMech,nodes);
+        dofrow = getLocalDoF(dofm,obj.fldMech,nodes);
 
         % kronecker delta in tensor form
         kron = [1;1;1;0;0;0];
@@ -130,10 +134,10 @@ function computeMatBiot(obj)
             Np = reshape(Nref',1,elem.nNode,nG);
             kron = repmat(kron,1,1,nG);
             iN = pagemtimes(kron,Np);
-            dofcol = getLocalDoF(obj.dofm,obj.fldFlow,nodes);
+            dofcol = getLocalDoF(dofm,obj.fldFlow,nodes);
           case "FVTPFA"
             iN = repmat(kron,1,1,nG);
-            dofcol = getLocalDoF(obj.dofm,obj.fldFlow,el);
+            dofcol = getLocalDoF(dofm,obj.fldFlow,el);
         end
 
         % compute local coupling matrix
@@ -163,8 +167,8 @@ function computeMatBiot(obj)
       uOld = getStateOld(obj,"displacements");
 
       % select active coefficients of solution vectors
-      entsPoro = obj.dofm.getActiveEntities(obj.fldMech,1);
-      entsFlow = obj.dofm.getActiveEntities(obj.fldFlow,1);
+      entsPoro = obj.domain.dofm.getActiveEntities(obj.fldMech,1);
+      entsFlow = obj.domain.dofm.getActiveEntities(obj.fldFlow,1);
 
       % get coupling blocks
       Qmech = getJacobian(obj,obj.fldMech,obj.fldFlow);

Code/discretizer/PhysicsSolvers/Flow/SinglePhaseFlow.m

@@ -35,10 +35,11 @@ function assembleSystem(obj,dt)
     end
 
     function updateState(obj,dSol)
+      dofm = obj.domain.dofm;
       if nargin > 1
-        ents = obj.dofm.getActiveEntities(obj.fieldId);
+        ents = dofm.getActiveEntities(obj.fieldId);
         state = getState(obj);
-        state.data.pressure(ents) = state.data.pressure(ents) + dSol(obj.dofm.getDoF(obj.fieldId));
+        state.data.pressure(ents) = state.data.pressure(ents) + dSol(dofm.getDoF(obj.fieldId));
       end
     end
 
@@ -73,12 +74,14 @@ function writeMatFile(obj,fac,tID)
     end
 
     function alpha = getRockCompressibility(obj,el)
-      if ismember(obj.mesh.cellTag(el),getTargetRegions(obj.dofm,["pressure","displacements"]))
+      mat = obj.domain.materials;
+      targetRegions = getTargetRegions(obj.domain.dofm,["pressure","displacements"]);
+      if ismember(obj.mesh.cellTag(el),targetRegions)
         alpha = 0; %this term is not needed in a coupled formulation
       else
-        if isfield(obj.materials.getMaterial(obj.mesh.cellTag(el)),"ConstLaw")
+        if isfield(mat.getMaterial(obj.mesh.cellTag(el)),"ConstLaw")
           %solid skeleton contribution to storage term as oedometric compressibility .
-          alpha = obj.materials.getMaterial(obj.mesh.cellTag(el)).ConstLaw.getRockCompressibility();
+          alpha = mat.getMaterial(obj.mesh.cellTag(el)).ConstLaw.getRockCompressibility();
         else
           alpha = 0;
         end
@@ -89,7 +92,7 @@ function writeMatFile(obj,fac,tID)
       % printPermeab - print the permeability for the cell or element.
       perm = zeros(obj.mesh.nCells,6);
       for el=1:obj.mesh.nCells
-        ktmp = obj.materials.getMaterial(obj.mesh.cellTag(el)).PorousRock.getPermMatrix();
+        ktmp = obj.domain.materials.getMaterial(obj.mesh.cellTag(el)).PorousRock.getPermMatrix();
         perm(el,1)=ktmp(1,1);
         perm(el,2)=ktmp(2,2);
         perm(el,3)=ktmp(3,3);