Add new save utility (+test case) for saving mortar projection stuff …

…without actually enforcing continuity.

fem/src/modules/SaveData/SaveProjection.F90

@@ -0,0 +1,195 @@
+!/*****************************************************************************/
+! *
+! *  Elmer, A Finite Element Software for Multiphysical Problems
+! *
+! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
+! * 
+! *  This program is free software; you can redistribute it and/or
+! *  modify it under the terms of the GNU General Public License
+! *  as published by the Free Software Foundation; either version 2
+! *  of the License, or (at your option) any later version.
+! * 
+! *  This program is distributed in the hope that it will be useful,
+! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! *  GNU General Public License for more details.
+! *
+! *  You should have received a copy of the GNU General Public License
+! *  along with this program (in file fem/GPL-2); if not, write to the 
+! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
+! *  Boston, MA 02110-1301, USA.
+! *
+! *****************************************************************************/
+!
+!/******************************************************************************
+! *
+! *  Subroutine for saving projections to new fields.
+! *
+! ******************************************************************************
+! *
+! *  Authors: Peter Råback
+! *  Email:   [email protected]
+! *  Web:     http://www.csc.fi/elmer
+! *  Address: CSC - IT Center for Science Ltd.
+! *           Keilaranta 14
+! *           02101 Espoo, Finland 
+! *
+! *  Original Date: 20 Nov 2001
+! *
+! *****************************************************************************/
+
+!> \ingroup Solvers
+!> \{
+
+
+!------------------------------------------------------------------------------
+SUBROUTINE SaveProjection_init( Model,Solver,dt,Transient )
+!------------------------------------------------------------------------------
+  USE DefUtils
+  IMPLICIT NONE
+!------------------------------------------------------------------------------
+  TYPE(Solver_t), TARGET :: Solver
+  TYPE(Model_t) :: Model
+  REAL(KIND=dp) :: dt
+  LOGICAL :: Transient
+!------------------------------------------------------------------------------  
+  CALL ListAddNewString( Solver % Values,'Variable',&
+      '-nooutput -global SaveProjection_var') 
+
+END SUBROUTINE SaveProjection_init
+!------------------------------------------------------------------------------
+
+
+!------------------------------------------------------------------------------
+!> Routine for saving projections from fields as fields. 
+!------------------------------------------------------------------------------
+SUBROUTINE SaveProjection( Model,Solver,dt,Transient )
+  USE DefUtils
+  USE Types
+  IMPLICIT NONE
+!------------------------------------------------------------------------------
+  TYPE(Solver_t), TARGET :: Solver
+  TYPE(Model_t) :: Model
+  REAL(KIND=dp) :: dt
+  LOGICAL :: Transient
+!------------------------------------------------------------------------------
+! Local variables
+!------------------------------------------------------------------------------
+  TYPE(ValueList_t), POINTER :: Params
+  TYPE(Variable_t), POINTER :: Var, TargetVar
+  INTEGER :: i,j,n,NoVar
+  CHARACTER(LEN=MAX_NAME_LEN) ::  VarName, TargetName
+  INTEGER, POINTER :: UnitPerm(:)
+  REAL(KIND=dp) :: Nrm
+  LOGICAL :: Found 
+
+  CALL Info('SaveProjection','Creating selected projected values as fields')
+
+  Params => GetSolverParams()
+
+  i = 0
+  DO WHILE(.TRUE.)  
+    i = i + 1
+    VarName = ListGetString( Params,'Variable '//I2S(i), Found )
+    IF(.NOT. Found) EXIT    
+  END DO
+  NoVar = i-1
+  CALL Info('SaveProjection','Saving projections from '//I2S(NoVar)//' fields')
+
+
+  DO i=1,NoVar
+    VarName = ListGetString( Params,'Variable '//I2S(i), Found )
+    Var => VariableGet( Model % Variables, TRIM(VarName) )
+    IF(.NOT. ASSOCIATED(Var)) THEN
+      CALL Warn('SaveProjection','Requested variable does not exist!')
+      CYCLE
+    END IF
+
+    TargetName = ListGetString( Params,'Target Variable '//I2S(i), Found )
+    IF(.NOT. Found) VarName = 'Projection '//TRIM(VarName)
+
+    TargetVar => VariableGet( Model % Variables, TRIM(TargetName) )
+    IF(.NOT. ASSOCIATED(TargetVar)) THEN
+      IF(.NOT. ASSOCIATED(Var % Perm) ) THEN
+        UnitPerm => NULL()
+        ALLOCATE(UnitPerm(SIZE(Var % Values)))
+        DO j=1,SIZE(Var % Values)
+          UnitPerm(j) = j
+        END DO
+        CALL VariableAddVector( Model % Mesh % Variables, Solver % Mesh, Solver, &
+            TRIM(TargetName), Var % Dofs, Perm = UnitPerm ) 
+      ELSE
+        CALL VariableAddVector( Model % Mesh % Variables, Solver % Mesh, Solver, &
+            TRIM(TargetName), Var % Dofs, Perm = Var % Perm, Secondary = .TRUE.)
+      END IF
+      TargetVar => VariableGet( Model % Variables, TRIM(TargetName) )       
+    END IF
+
+    ! Do additive projection!
+    CALL ProjectToVariable()
+    Nrm = Nrm + SUM(TargetVar % Values**2)
+  END DO
+
+  Nrm = SQRT(Nrm)
+  IF(SIZE(Solver % Variable % Values) == 1 ) THEN
+    Solver % Variable % Values = Nrm
+  END IF
+
+  WRITE(Message,'(A,ES12.3)') 'Combined L2 norm of all projected fields: ', Nrm
+  CALL Info('SaveProjection',Message)
+
+CONTAINS
+
+
+  SUBROUTINE ProjectToVariable()
+    TYPE(Matrix_t), POINTER :: A
+    INTEGER :: bc, dofs, i, j, k, pi, pj
+    INTEGER, POINTER :: Rows(:), Cols(:)
+    REAL(KIND=dp) :: r1
+    REAL(KIND=dp), POINTER :: Values(:)      
+
+    dofs = Var % Dofs
+    TargetVar % Values = 0.0_dp
+
+    ! Go through all the projectors.
+    ! There could be perhaps reason to skip some, but this will do for now.
+    DO bc=1,Model % NumberOfBCs        
+      A => CurrentModel % BCs(bc) % PMatrix
+      IF(.NOT. ASSOCIATED(A) ) THEN
+        A => Solver % MortarBCs(bc) % Projector
+      END IF
+
+      IF(.NOT. ASSOCIATED(A)) CYCLE
+      n = A % NumberOfRows
+      Rows => A % Rows
+      Cols => A % Cols
+      Values => A % Values
+
+      DO k=1,dofs
+        DO i=1,n
+          IF(ASSOCIATED(A % InvPerm)) THEN
+            pi = A % InvPerm(i)
+          ELSE
+            pi = i
+          END IF
+          IF(ASSOCIATED(Var % Perm)) pi = Var % Perm(pi)          
+          IF(pi==0) CYCLE
+          pi = dofs*(pi-1)+k
+          r1 = 0.0_dp
+          DO j=Rows(i),Rows(i+1)-1
+            pj = Cols(j)
+            IF(ASSOCIATED(Var % Perm)) pj = Var % Perm(pj)
+            IF(pj==0) CYCLE
+            pj = dofs*(pj-1)+k
+            r1 = r1 + Values(j) * Var % Values(pj)
+          END DO
+
+          TargetVar % Values(pi) = TargetVar % Values(pi) + r1 
+        END DO
+      END DO
+    END DO
+  END SUBROUTINE ProjectToVariable
+
+END SUBROUTINE SaveProjection
+
+!> \}

fem/tests/MortarSaveProjection/CMakeLists.txt

@@ -0,0 +1,8 @@
+INCLUDE(test_macros)
+INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}/fem/src)
+
+CONFIGURE_FILE( case.sif case.sif COPYONLY)
+
+file(COPY squares.grd ELMERSOLVER_STARTINFO DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
+
+ADD_ELMER_TEST(MortarSaveProjection LABELS quick mortar)

fem/tests/MortarSaveProjection/ELMERSOLVER_STARTINFO

@@ -0,0 +1 @@
+case.sif

fem/tests/MortarSaveProjection/case.sif

@@ -0,0 +1,137 @@
+! Test case for just applying projectors to map stuff without enforcing
+! any kind of continuity.
+!
+! P.R. 1.10.2024
+
+
+Header
+  CHECK KEYWORDS Warn
+  Mesh DB "." "squares"
+  Include Path ""
+  Results Directory ""
+End
+
+Simulation
+  Max Output Level = 5
+  Coordinate System = "Cartesian 2D"
+  Coordinate Mapping(3) = 1 2 3
+  Simulation Type = Steady State
+  Steady State Max Iterations = 1
+  Output Intervals = 1
+
+  Output Intervals = 0
+
+!  Post File = case.vtu
+!  Ascii Output = True
+End 
+
+Constants
+  Gravity(4) = 0 -1 0 9.82
+  Stefan Boltzmann = 5.67e-08
+  Permittivity of Vacuum = 8.8542e-12
+  Boltzmann Constant = 1.3807e-23
+  Unit Charge = 1.602e-19
+End
+
+Body 1
+  Target Bodies(1) = 1
+  Name = "Body1"
+  Equation = 1
+  Material = 1
+  Body Force = 1
+End
+
+Body 2
+  Target Bodies(1) = 2
+  Name = "Body2"
+  Equation = 1
+  Material = 1
+End
+
+Solver 1
+  Equation = Heat Equation
+  Procedure = "HeatSolve" "HeatSolver"
+  Variable = Temperature
+
+  Nonlinear System Max Iterations = 1
+
+  Linear System Solver = Iterative
+  Linear System Iterative Method = BiCGStabl
+  Linear System Max Iterations = 5000
+  Linear System Convergence Tolerance = 1.0e-8
+  Linear System Preconditioning = ILU1
+  Linear System Residual Output = 20
+  Linear System Precondition Recompute = 1
+
+! We do not apply, only do afterwords projection.
+!  Apply Mortar BCs = Logical True
+End
+
+
+Solver 2
+! Settings mainly for timing and verification
+!  Exec Solver = never
+
+  Equation = SaveProject
+  Procedure = "SaveData" "SaveProjection"
+
+  Variable 1 = Temperature
+  Variable 2 = Coordinate 1
+  Variable 3 = Coordinate 2
+
+! Note that the field values depend on type of projection.
+! They may, or may not, be multiplied by nodal weight. 
+  Target Variable 1 = Proj Temperature
+  Target Variable 2 = Proj X
+  Target Variable 3 = Proj Y 
+
+  Apply Mortar BCs = Logical True
+End 
+
+
+Equation 1
+  Name = "Heat"
+  Active Solvers(2) = 1 2
+End
+
+Material 1
+  Name = "Ideal"
+  Heat Conductivity = 1
+  Heat Capacity = 1
+  Density = 1
+End
+
+Body Force 1
+  Name = "Heating"
+  Heat Source = 1.0
+End
+
+Boundary Condition 1
+  Target Boundaries(1) = 4
+  Name = "Left-left"
+  Temperature = 1.0
+End
+
+Boundary Condition 2
+  Target Boundaries(1) = 6
+  Name = "Right-Right"
+  Temperature = 2.0
+End
+
+Boundary Condition 3
+  Target Boundaries(1) = 2
+  Name = "Left-Right"
+  Mortar BC = Integer 4 
+  Galerkin Projector = Logical True
+  Plane Projector = Logical True
+End
+
+Boundary Condition 4
+  Target Boundaries(1) = 8
+  Name = "Right-Left"  
+End
+
+
+Solver 1 :: Reference Norm = 1.81336882E+00
+Solver 2 :: Reference Norm = 1.57162326E-01
+

fem/tests/MortarSaveProjection/runtest.cmake

@@ -0,0 +1,3 @@
+include(test_macros)
+execute_process(COMMAND ${ELMERGRID_BIN} 1 2 squares.grd)
+RUN_ELMER_TEST()

fem/tests/MortarSaveProjection/squares.grd

@@ -0,0 +1,31 @@
+***** ElmerGrid input file for structured grid generation *****
+Version = 210903
+Coordinate System = Cartesian 2D
+Subcell Divisions in 2D = 3 2
+Subcell Sizes 1 = 1 0.1 1.5 
+Subcell Sizes 2 = 1 0.5 
+Material Structure in 2D
+  0 0 2
+  1 0 2 
+End
+
+Element Degree = 1
+Element Innernodes = False
+
+Materials Interval = 1 2
+Boundary Definitions
+! type     out      int
+   1      -1        1        1
+   2      -2        1        1
+   3      -3        1        1
+   4      -4        1        1
+   5      -1        2        1
+   6      -2        2        1
+   7      -3        2        1
+   8      -4        2        1
+End
+Element Ratios 1 = 1 1 1 
+Element Ratios 2 = 1 1 
+Element Divisions 1 = 10 0 15 
+Element Divisions 2 = 10 5 
+