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 
+