Re-work application of Dirichlet boundary conditions to work with multiple...

Re-work application of Dirichlet boundary conditions to work with multiple time function for the same dof.

+mafe/FeProblem.m

@@ -150,21 +150,19 @@ classdef FeProblem < handle
         end
         %% apply constraints of type Neumann to the rhs
         function [b, on, off] = applyConstraintsNeumann(self, b, tfun)
+            % Default value for tfun (if not given) is mafe.TimeFunction.Static
+            if ~exist('tfun', 'var')
+              tfun = mafe.TimeFunction.Static;
+            end
             % apply constraints - dynamic/static
             idx = [];
             val = [];
             for con = self.const
                 if con.type == mafe.ConstraintType.Neumann
-                    idx = [ idx, con.getDofSysIndices() ];
-                    if ~exist('tfun', 'var')
-                        % tfun not given -> time-constant dynamic/static constraint
-                        val = [ val, con.value() ];
-                    elseif con.tfun == tfun
+                    if con.tfun == tfun
                         % collect matching tfun value
+                        idx = [ idx, con.getDofSysIndices() ];
                         val = [ val, con.value() ];
-                    else
-                        % collect zero (vector value)
-                        val = [ val, 0.0*con.value() ];
                     end
                 end
             end
@@ -177,25 +175,41 @@ classdef FeProblem < handle
             % create index vector of dof system indices with    Neumann constraints
             on  = idx;
         end
-        %% apply constraints of type Dirichlet to coefficent matrix and rhs
+        %% apply constraints of type Dirichlet to coefficient matrix and rhs
         function [b, on, off, A] = applyConstraintsDirichlet(self, A, b, tfun)
+            % Default value for tfun (if not given) is mafe.TimeFunction.Static
+            if ~exist('tfun', 'var')
+              tfun = mafe.TimeFunction.Static;
+            end
             % apply constraints - kinematic
+            % --- pass 1: (homogeneous case)
+            idx_h = [];
+            for con = self.const
+              if con.type == mafe.ConstraintType.Dirichlet
+                % constraint of a different tfun is treated as homogeneous
+                if con.tfun ~= tfun
+                  idx_h = [ idx_h, con.getDofSysIndices() ];
+                end
+              end
+            end
+            % check for repeated indices in kinematic constraints
+            if length(idx_h) ~= length(unique(idx_h))
+               %disp('WARNING: FeProblem::applyConstraints() : Repeated indices in Dirichlet constraint !');
+               % give preference to first occurence
+               [~,occurence,~] = unique(idx_h, 'first');
+               idx_h = idx_h(occurence);
+            end
+            % --- pass 2: (inhomogeneous case)
             idx = [];
             val = [];
             for con = self.const
-                if con.type == mafe.ConstraintType.Dirichlet
-                    idx = [ idx, con.getDofSysIndices() ];
-                    if ~exist('tfun', 'var')
-                        % tfun not given -> time-constant kinematic constraint
-                        val = [ val, con.value() ];
-                    elseif con.tfun == tfun
-                        % collect matching tfun value
-                        val = [ val, con.value() ];
-                    else
-                        % collect zero, treat this as homogeneous constraint
-                        val = [ val, 0.0*con.value() ];
-                    end
+              if con.type == mafe.ConstraintType.Dirichlet
+                % constraint matches given tfun
+                if con.tfun == tfun
+                  idx = [ idx, con.getDofSysIndices() ];
+                  val = [ val, con.value() ];
                 end
+              end
             end
             % check for repeated indices in kinematic constraints
             if length(idx) ~= length(unique(idx))
@@ -206,9 +220,16 @@ classdef FeProblem < handle
                idx = idx(occurence);
             end
             % impose inhomogeneous Dirchlet conditions
-            b = b - A(:,idx) * val'; % assuming distinct indices here
+            if length(idx)>0
+              b = b - A(:,idx) * val'; % assuming distinct indices here
+            end
+            % put rhs values
+            b(idx_h) = 0.0; % homogeneous first
+            b(idx  ) = val; % inhomogeneous second
+            % determine affected indices
+            idx = unique([idx_h, idx]);
+            % modify lhs
             A(idx,:) = 0.0; A(:,idx) = 0.0; A(idx,idx) = eye(length(idx));
-            b(idx) = val;
             % create index vector of dof system indices without Dirichlet constraints
             off = setdiff([1:self.ndofs], idx);
             % create index vector of dof system indices with    Dirichlet constraints