diff --git a/pages/BaseMethods.md b/pages/BaseMethods.md
index b1c7c6c09..5b586f4e9 100644
--- a/pages/BaseMethods.md
+++ b/pages/BaseMethods.md
@@ -1,6 +1,6 @@
 # BaseMethods
 
-`BaseMethods` library contains the modules which defines and implements methods  (routines) for data types defined in BaseType.
+`BaseMethods` library contains the modules which defines and implements methods (routines) for data types defined in BaseType.
 
 At present BaseMethods contains following modules.
 
diff --git a/src/modules/BaseInterpolation/src/BaseInterpolation_Method.F90 b/src/modules/BaseInterpolation/src/BaseInterpolation_Method.F90
index 4afe02b7f..c7470dccc 100644
--- a/src/modules/BaseInterpolation/src/BaseInterpolation_Method.F90
+++ b/src/modules/BaseInterpolation/src/BaseInterpolation_Method.F90
@@ -22,17 +22,15 @@ MODULE BaseInterpolation_Method
 USE String_Class, ONLY: String
 USE StringUtility, ONLY: UpperCase
 USE Display_Method, ONLY: Tostring
-USE BaseType, ONLY: poly => TypePolynomialOpt, &
-                    ip => TypeQuadratureOpt, &
-                    BaseInterpolation_, &
-                    LagrangeInterpolation_, &
-                    SerendipityInterpolation_, &
-                    HermitInterpolation_, &
-                    HierarchyInterpolation_, &
-                    OrthogonalInterpolation_
-
+USE BaseType, ONLY: poly => TypePolynomialOpt
+USE BaseType, ONLY: ip => TypeQuadratureOpt
+USE BaseType, ONLY: BaseInterpolation_
+USE BaseType, ONLY: LagrangeInterpolation_
+USE BaseType, ONLY: SerendipityInterpolation_
+USE BaseType, ONLY: HermitInterpolation_
+USE BaseType, ONLY: HierarchyInterpolation_
+USE BaseType, ONLY: OrthogonalInterpolation_
 IMPLICIT NONE
-
 PRIVATE
 
 PUBLIC :: ASSIGNMENT(=)
@@ -226,8 +224,11 @@ FUNCTION InterpolationPoint_ToInteger(name) RESULT(ans)
 
   SELECT CASE (astr)
 
+  CASE ("CENTER")
+    ans = ip%Center
+
   CASE ("EQUIDISTANCE")
-    ans = ip%equidistance
+    ans = ip%Equidistance
 
   CASE ("GAUSSLEGENDRE")
     ans = ip%GaussLegendre
@@ -580,6 +581,13 @@ FUNCTION InterpolationPoint_ToChar(name, isUpper) RESULT(ans)
   IF (PRESENT(isUpper)) isUpper0 = isUpper
 
   SELECT CASE (name)
+  CASE (ip%Center)
+    IF (isUpper0) THEN
+      ans = "CENTER"
+    ELSE
+      ans = "Center"
+    END IF
+
   CASE (ip%equidistance)
     IF (isUpper0) THEN
       ans = "EQUIDISTANCE"
diff --git a/src/modules/BaseMethod/src/BaseMethod.F90 b/src/modules/BaseMethod/src/BaseMethod.F90
index 8f30d9e34..23361c401 100644
--- a/src/modules/BaseMethod/src/BaseMethod.F90
+++ b/src/modules/BaseMethod/src/BaseMethod.F90
@@ -99,7 +99,6 @@ MODULE BaseMethod
 USE IntVector_Method
 USE IndexValue_Method
 USE KeyValue_Method
-USE IterationData_Method
 USE Vector3D_Method
 USE RealVector_Method
 USE DOF_Method
diff --git a/src/modules/BaseType/src/BaseType.F90 b/src/modules/BaseType/src/BaseType.F90
index a7ee2020e..6b8c5d1fb 100644
--- a/src/modules/BaseType/src/BaseType.F90
+++ b/src/modules/BaseType/src/BaseType.F90
@@ -33,10 +33,10 @@ MODULE BaseType
 
 USE GlobalData, ONLY: RelativeConvergence, ConvergenceInRes, &
                       ConvergenceInSol, ConvergenceInResSol, &
-                      AbsoluteConvergence, NormL2, &
+                      AbsoluteConvergence, NormL1, NormL2, NormInfinity, &
                       StressTypeVoigt, OMP_THREADS_JOINED
 
-USE GlobalData, ONLY: Equidistance, EquidistanceQP, GaussQP, &
+USE GlobalData, ONLY: Equidistance, CenterQP, EquidistanceQP, GaussQP, &
                       GaussLegendreQP, GaussLegendreLobattoQP, &
                       GaussLegendreRadau, GaussLegendreRadauLeft, &
                       GaussLegendreRadauRight, GaussRadauQP, &
@@ -141,9 +141,6 @@ MODULE BaseType
 PUBLIC :: CSRMatrix_
 PUBLIC :: TypeCSRMatrix
 PUBLIC :: CSRMatrixPointer_
-PUBLIC :: IterationData_
-PUBLIC :: TypeIterationData
-PUBLIC :: IterationDataPointer_
 PUBLIC :: VoigtRank2Tensor_
 PUBLIC :: TypeVoigtRank2Tensor
 PUBLIC :: VoigtRank2TensorPointer
@@ -312,6 +309,8 @@ MODULE BaseType
   REAL(DFP) :: minus_one = -1.0_DFP
   REAL(DFP) :: two = 2.0_DFP
   REAL(DFP) :: minus_two = -2.0_DFP
+  REAL(DFP) :: three = 3.0_DFP
+  REAL(DFP) :: minus_three = -3.0_DFP
   REAL(DFP) :: pi = 3.14159265359_DFP
   REAL(DFP) :: two_pi = 2.0_DFP * 3.14159265359_DFP
   REAL(DFP) :: pi_by_two = 0.5_DFP * 3.14159265359_DFP
@@ -439,7 +438,7 @@ MODULE BaseType
 
 TYPE :: RealVector_
   INTEGER(I4B) :: tDimension = 1_I4B
-  REAL(DFP), ALLOCATABLE :: Val(:)
+  REAL(DFP), ALLOCATABLE :: val(:)
 END TYPE RealVector_
 
 TYPE(RealVector_), PARAMETER :: TypeRealVector = RealVector_(Val=NULL())
@@ -534,6 +533,13 @@ MODULE BaseType
 ! summary: Degree of freedom object type
 
 TYPE :: DOF_
+  INTEGER(I4B) :: storageFMT = FMT_NODES
+    !! Storage format
+  INTEGER(I4B) :: mapRow = 0
+    !! Number of rows in map which contains useful data
+    !! Number of physical variables = mapRow - 1
+  INTEGER(I4B) :: valMapSize = 0
+    !! The size of valMap which contains useful data
   INTEGER(I4B), ALLOCATABLE :: map(:, :)
     !! Encapsulation of information of DOF
     !! map contains 6 columns
@@ -550,13 +556,6 @@ MODULE BaseType
     !!    For example, map(n+1, 4) contains the total DOF
   INTEGER(I4B), ALLOCATABLE :: valMap(:)
     !! Val map
-  INTEGER(I4B) :: storageFMT = FMT_NODES
-    !! Storage format
-  INTEGER(I4B) :: mapRow = 0
-    !! Number of rows in map which contains useful data
-    !! Number of physical variables = mapRow - 1
-  INTEGER(I4B) :: valMapSize = 0
-    !! The size of valMap which contains useful data
 END TYPE DOF_
 
 TYPE(DOF_), PARAMETER :: TypeDOF = DOF_(map=NULL(), valMap=NULL())
@@ -705,65 +704,6 @@ MODULE BaseType
   CLASS(CSRMatrix_), POINTER :: ptr => NULL()
 END TYPE CSRMatrixPointer_
 
-!----------------------------------------------------------------------------
-!                                                            IterationData_
-!----------------------------------------------------------------------------
-
-!> author: Vikas Sharma, Ph. D.
-! date: 14 June 2022
-! summary: Iteration data
-
-TYPE :: IterationData_
-  INTEGER(I4B) :: maxIter = 100
-    !! Maximum number of iterations allowed
-  INTEGER(I4B) :: iterationNumber = 1
-    !! Iteration number
-  REAL(DFP) :: residualError0 = 0.0
-    !! Initial Residual error
-  REAL(DFP) :: residualError = 0.0
-    !! Current residual error
-  REAL(DFP) :: residualTolerance = 1.0E-5
-    !! Tolerance for checking convergence in residual
-  REAL(DFP) :: solutionError0 = 0.0
-    !! Initial solution error
-  REAL(DFP) :: solutionError = 0.0
-    !! Current solution error
-  REAL(DFP) :: solutionTolerance = 1.0E-5
-    !! Tolerance for checking convergence in solution
-  INTEGER(I4B) :: convergenceType = RelativeConvergence
-    !! Type of convergence
-  INTEGER(I4B) :: convergenceIn = ConvergenceInRes
-    !! Check Convergence in solution and/or residual
-  INTEGER(I4B) :: normType = NormL2
-    !! Error norm type
-  LOGICAL(LGT) :: converged = .FALSE.
-    !! Status of convergence
-  REAL(DFP) :: timeAtStart = 0.0
-    !! Starting time
-  REAL(DFP) :: timeAtEnd = 0.0
-    !! Present time
-  REAL(DFP), ALLOCATABLE :: convergenceData(:, :)
-    !! history of convergence data
-    !! each column corresponding to a iteration
-  TYPE(String), ALLOCATABLE :: header(:)
-    !! header for convergenceData
-END TYPE IterationData_
-
-!----------------------------------------------------------------------------
-!
-!----------------------------------------------------------------------------
-
-TYPE(IterationData_), PARAMETER :: TypeIterationData = &
-                                   IterationData_(header=NULL())
-
-!----------------------------------------------------------------------------
-!
-!----------------------------------------------------------------------------
-
-TYPE :: IterationDataPointer_
-  CLASS(IterationData_), POINTER :: ptr => NULL()
-END TYPE IterationDataPointer_
-
 !----------------------------------------------------------------------------
 !                                                       VoigtRank2Tensor_
 !----------------------------------------------------------------------------
@@ -2171,6 +2111,9 @@ END SUBROUTINE InterfaceMatrixSubroutine
   INTEGER(I4B) :: both = convergenceInResSol
   INTEGER(I4B) :: relative = relativeConvergence
   INTEGER(I4B) :: absolute = absoluteConvergence
+  INTEGER(I4B) :: normL1 = NormL1
+  INTEGER(I4B) :: normL2 = NormL2
+  INTEGER(I4B) :: normInfinity = NormInfinity
 END TYPE ConvergenceOpt_
 
 TYPE(ConvergenceOpt_), PARAMETER :: TypeConvergenceOpt = ConvergenceOpt_()
@@ -2472,6 +2415,7 @@ END SUBROUTINE InterfaceMatrixSubroutine
   INTEGER(I4B) :: BlythPozChebyshev = BlythPozChebyshevQP
   INTEGER(I4B) :: IsaacLegendre = IsaacLegendreQP
   INTEGER(I4B) :: IsaacChebyshev = IsaacChebyshevQP
+  INTEGER(I4B) :: Center = CenterQP
   INTEGER(I4B) :: default = GaussLegendreQP
 END TYPE QuadratureOpt_
 
diff --git a/src/modules/CMakeLists.txt b/src/modules/CMakeLists.txt
index 57ea10c0d..6ec99ca74 100644
--- a/src/modules/CMakeLists.txt
+++ b/src/modules/CMakeLists.txt
@@ -140,9 +140,6 @@ include(${CMAKE_CURRENT_LIST_DIR}/IntVector/CMakeLists.txt)
 # IndexValue
 include(${CMAKE_CURRENT_LIST_DIR}/IndexValue/CMakeLists.txt)
 
-# IndexValue
-include(${CMAKE_CURRENT_LIST_DIR}/IterationData/CMakeLists.txt)
-
 # KeyValue
 include(${CMAKE_CURRENT_LIST_DIR}/KeyValue/CMakeLists.txt)
 
@@ -227,6 +224,9 @@ include(${CMAKE_CURRENT_LIST_DIR}/CSRSparsity/CMakeLists.txt)
 # CSRMatrix
 include(${CMAKE_CURRENT_LIST_DIR}/CSRMatrix/CMakeLists.txt)
 
+# ConvergenceOpt
+include(${CMAKE_CURRENT_LIST_DIR}/ConvergenceOpt/CMakeLists.txt)
+
 # BaseMethod
 include(${CMAKE_CURRENT_LIST_DIR}/BaseMethod/CMakeLists.txt)
 
diff --git a/src/modules/ConvergenceOpt/CMakeLists.txt b/src/modules/ConvergenceOpt/CMakeLists.txt
new file mode 100644
index 000000000..26f45560d
--- /dev/null
+++ b/src/modules/ConvergenceOpt/CMakeLists.txt
@@ -0,0 +1,19 @@
+# This program is a part of EASIFEM library
+# Copyright (C) 2020-2021  Vikas Sharma, Ph.D
+#
+# 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 3 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.  If not, see <https: //www.gnu.org/licenses/>
+#
+
+set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/")
+target_sources(${PROJECT_NAME} PRIVATE ${src_path}/ConvergenceOptUtility.F90)
diff --git a/src/modules/ConvergenceOpt/src/ConvergenceOptUtility.F90 b/src/modules/ConvergenceOpt/src/ConvergenceOptUtility.F90
new file mode 100644
index 000000000..b91af6270
--- /dev/null
+++ b/src/modules/ConvergenceOpt/src/ConvergenceOptUtility.F90
@@ -0,0 +1,180 @@
+! This program is a part of EASIFEM library
+! Expandable And Scalable Infrastructure for Finite Element Methods
+! htttps://www.easifem.com
+!
+! 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 3 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.  If not, see <https: //www.gnu.org/licenses/>
+!
+
+MODULE ConvergenceOptUtility
+USE GlobalData, ONLY: DFP, I4B, LGT
+USE String_Class, ONLY: String
+IMPLICIT NONE
+PRIVATE
+
+PUBLIC :: ConvergenceType_ToInt
+PUBLIC :: ConvergenceType_ToChar
+PUBLIC :: ConvergenceType_ToString
+PUBLIC :: ConvergenceIn_ToInt
+PUBLIC :: ConvergenceIn_ToChar
+PUBLIC :: ConvergenceIn_ToString
+PUBLIC :: NormType_ToInt
+PUBLIC :: NormType_ToChar
+PUBLIC :: NormType_ToString
+
+!----------------------------------------------------------------------------
+!                                                       ConvergenceType_ToInt
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from character to integer
+
+INTERFACE
+  MODULE FUNCTION ConvergenceType_ToInt(name) RESULT(ans)
+    CHARACTER(*), INTENT(IN) :: name
+    INTEGER(I4B) :: ans
+  END FUNCTION ConvergenceType_ToInt
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                    ConvergenceType_ToChar
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from integer to character
+
+INTERFACE
+  MODULE FUNCTION ConvergenceType_ToChar(name, isUpper) RESULT(ans)
+    INTEGER(I4B), INTENT(IN) :: name
+    LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isUpper
+    CHARACTER(:), ALLOCATABLE :: ans
+  END FUNCTION ConvergenceType_ToChar
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                    ConvergenceType_ToString
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from integer to character
+
+INTERFACE
+  MODULE FUNCTION ConvergenceType_ToString(name, isUpper) RESULT(ans)
+    INTEGER(I4B), INTENT(IN) :: name
+    LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isUpper
+    TYPE(String) :: ans
+  END FUNCTION ConvergenceType_ToString
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                        ConvergenceIn_ToInt
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from character to integer
+
+INTERFACE
+  MODULE FUNCTION ConvergenceIn_ToInt(name) RESULT(ans)
+    CHARACTER(*), INTENT(IN) :: name
+    INTEGER(I4B) :: ans
+  END FUNCTION ConvergenceIn_ToInt
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                       ConvergenceIn_ToChar
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from integer to character
+
+INTERFACE
+  MODULE FUNCTION ConvergenceIn_ToChar(name, isUpper) RESULT(ans)
+    INTEGER(I4B), INTENT(IN) :: name
+    LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isUpper
+    CHARACTER(:), ALLOCATABLE :: ans
+  END FUNCTION ConvergenceIn_ToChar
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                     ConvergenceIn_ToString
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from integer to character
+
+INTERFACE
+  MODULE FUNCTION ConvergenceIn_ToString(name, isUpper) RESULT(ans)
+    INTEGER(I4B), INTENT(IN) :: name
+    LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isUpper
+    TYPE(String) :: ans
+  END FUNCTION ConvergenceIn_ToString
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                             NormType_ToInt
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from character to integer
+
+INTERFACE
+  MODULE FUNCTION NormType_ToInt(name) RESULT(ans)
+    CHARACTER(*), INTENT(IN) :: name
+    INTEGER(I4B) :: ans
+  END FUNCTION NormType_ToInt
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                            NormType_ToChar
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from integer to character
+
+INTERFACE
+  MODULE FUNCTION NormType_ToChar(name, isUpper) RESULT(ans)
+    INTEGER(I4B), INTENT(IN) :: name
+    LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isUpper
+    CHARACTER(:), ALLOCATABLE :: ans
+  END FUNCTION NormType_ToChar
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!                                                           NormType_ToString
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2026-07-14
+! summary: convert convergence from integer to character
+
+INTERFACE
+  MODULE FUNCTION NormType_ToString(name, isUpper) RESULT(ans)
+    INTEGER(I4B), INTENT(IN) :: name
+    LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isUpper
+    TYPE(String) :: ans
+  END FUNCTION NormType_ToString
+END INTERFACE
+
+!----------------------------------------------------------------------------
+!
+!----------------------------------------------------------------------------
+
+END MODULE ConvergenceOptUtility
diff --git a/src/modules/ElemshapeData/src/ElemshapeData_InterpolMethods.F90 b/src/modules/ElemshapeData/src/ElemshapeData_InterpolMethods.F90
index b76509037..926ed7497 100644
--- a/src/modules/ElemshapeData/src/ElemshapeData_InterpolMethods.F90
+++ b/src/modules/ElemshapeData/src/ElemshapeData_InterpolMethods.F90
@@ -15,11 +15,13 @@
 ! along with this program.  If not, see <https: //www.gnu.org/licenses/>
 !
 !
-! This file contains the interpolation methods interfaces\
+! This file contains the interpolation methods interfaces
 
 MODULE ElemshapeData_InterpolMethods
 USE GlobalData, ONLY: DFP, I4B, LGT
-USE BaseType, ONLY: ElemShapeData_, STElemShapeData_, FEVariable_
+USE BaseType, ONLY: ElemShapeData_
+USE BaseType, ONLY: STElemShapeData_
+USE BaseType, ONLY: FEVariable_
 IMPLICIT NONE
 PRIVATE
 
@@ -35,7 +37,7 @@ MODULE ElemshapeData_InterpolMethods
 ! date: 4 March 2021
 ! summary: returns the interpolation of a FEVariable
 !
-!# Introduction
+!# GetInterpolation
 !
 ! If ans is not initiated then it will be initiated
 ! If ans is initiated then we will just call GetInterpolation_
@@ -59,16 +61,12 @@ MODULE ElemshapeData_InterpolMethods
 !                                we do not have time shape function for
 !                                all quadrature points in time in obj)
 
-INTERFACE
+INTERFACE GetInterpolation
   MODULE PURE SUBROUTINE GetInterpolation1(obj, ans, val)
     CLASS(ElemshapeData_), INTENT(IN) :: obj
     TYPE(FEVariable_), INTENT(INOUT) :: ans
     TYPE(FEVariable_), INTENT(IN) :: val
   END SUBROUTINE GetInterpolation1
-END INTERFACE
-
-INTERFACE GetInterpolation
-  MODULE PROCEDURE GetInterpolation1
 END INTERFACE GetInterpolation
 
 !----------------------------------------------------------------------------
@@ -79,7 +77,7 @@ END SUBROUTINE GetInterpolation1
 ! date: 4 March 2021
 ! summary: returns the interpolation of a FEVariable
 !
-!# Introduction
+!# GetInterpolation_
 !
 ! - Returns the interpolation of a FEVariable_
 ! - The result is returned in ans
@@ -91,16 +89,12 @@ END SUBROUTINE GetInterpolation1
 ! - the val can be defined on quadrature (do nothing) or nodal (interpol)
 ! - The `vartype` of val can be constant, space, time, spacetime
 
-INTERFACE
+INTERFACE GetInterpolation_
   MODULE PURE SUBROUTINE GetInterpolation_1(obj, ans, val)
     CLASS(ElemshapeData_), INTENT(IN) :: obj
     TYPE(FEVariable_), INTENT(INOUT) :: ans
     TYPE(FEVariable_), INTENT(IN) :: val
   END SUBROUTINE GetInterpolation_1
-END INTERFACE
-
-INTERFACE GetInterpolation_
-  MODULE PROCEDURE GetInterpolation_1
 END INTERFACE GetInterpolation_
 
 !----------------------------------------------------------------------------
@@ -111,7 +105,7 @@ END SUBROUTINE GetInterpolation_1
 ! date: 4 March 2021
 ! summary: returns the interpolation of a FEVariable
 !
-!# Introduction
+!# GetInterpolation_
 !
 ! - Returns the interpolation of a FEVariable_
 ! - The result is returned in ans
@@ -123,7 +117,7 @@ END SUBROUTINE GetInterpolation_1
 ! - the val can be defined on quadrature (do nothing) or nodal (interpol)
 ! - The `vartype` of val can be constant, space, time, spacetime
 !
-INTERFACE
+INTERFACE GetInterpolation_
   MODULE PURE SUBROUTINE GetInterpolation_1a(obj, ans, val, scale, &
                                              addContribution)
     CLASS(ElemshapeData_), INTENT(IN) :: obj
@@ -132,10 +126,6 @@ MODULE PURE SUBROUTINE GetInterpolation_1a(obj, ans, val, scale, &
     REAL(DFP), INTENT(IN) :: scale
     LOGICAL, INTENT(IN) :: addContribution
   END SUBROUTINE GetInterpolation_1a
-END INTERFACE
-
-INTERFACE GetInterpolation_
-  MODULE PROCEDURE GetInterpolation_1a
 END INTERFACE GetInterpolation_
 
 !----------------------------------------------------------------------------
@@ -146,7 +136,7 @@ END SUBROUTINE GetInterpolation_1a
 ! date: 4 March 2021
 ! summary: returns the interpolation of a FEVariable
 !
-!# Introduction
+!# GetInterpolation
 !
 ! If ans is not initiated then it will be initiated. If
 ! ans is initiated then its properties will not be altered.
@@ -162,16 +152,12 @@ END SUBROUTINE GetInterpolation_1a
 !
 ! - ans will Quadrature and SpaceTime
 
-INTERFACE
+INTERFACE GetInterpolation
   MODULE PURE SUBROUTINE GetInterpolation2(obj, ans, val)
     CLASS(STElemshapeData_), INTENT(IN) :: obj(:)
     TYPE(FEVariable_), INTENT(INOUT) :: ans
     TYPE(FEVariable_), INTENT(IN) :: val
   END SUBROUTINE GetInterpolation2
-END INTERFACE
-
-INTERFACE GetInterpolation
-  MODULE PROCEDURE GetInterpolation2
 END INTERFACE GetInterpolation
 
 !----------------------------------------------------------------------------
@@ -193,17 +179,13 @@ END SUBROUTINE GetInterpolation2
 ! - The val can have following ranks; scalar, vector, matrix
 ! - the val can be defined on quadrature (do nothing) or nodal (interpol)
 ! - The `vartype` of val can be constant, space, time, spacetime
-!
-INTERFACE
+
+INTERFACE GetInterpolation_
   MODULE PURE SUBROUTINE GetInterpolation_2(obj, ans, val)
     CLASS(STElemshapeData_), INTENT(IN) :: obj(:)
     TYPE(FEVariable_), INTENT(INOUT) :: ans
     TYPE(FEVariable_), INTENT(IN) :: val
   END SUBROUTINE GetInterpolation_2
-END INTERFACE
-
-INTERFACE GetInterpolation_
-  MODULE PROCEDURE GetInterpolation_2
 END INTERFACE GetInterpolation_
 
 !----------------------------------------------------------------------------
@@ -226,7 +208,7 @@ END SUBROUTINE GetInterpolation_2
 ! - the val can be defined on quadrature (do nothing) or nodal (interpol)
 ! - The `vartype` of val can be constant, space, time, spacetime
 
-INTERFACE
+INTERFACE GetInterpolation_
   MODULE PURE SUBROUTINE GetInterpolation_2a(obj, ans, val, scale, &
                                              addContribution)
     CLASS(STElemshapeData_), INTENT(IN) :: obj(:)
@@ -235,10 +217,6 @@ MODULE PURE SUBROUTINE GetInterpolation_2a(obj, ans, val, scale, &
     REAL(DFP), INTENT(IN) :: scale
     LOGICAL, INTENT(IN) :: addContribution
   END SUBROUTINE GetInterpolation_2a
-END INTERFACE
-
-INTERFACE GetInterpolation_
-  MODULE PROCEDURE GetInterpolation_2a
 END INTERFACE GetInterpolation_
 
 !----------------------------------------------------------------------------
@@ -250,16 +228,16 @@ END SUBROUTINE GetInterpolation_2a
 ! update: 2021-12-13
 ! summary: Interpolation of FEVariable
 
-INTERFACE
+INTERFACE Interpolation
   MODULE PURE FUNCTION Interpolation1(obj, val) RESULT(ans)
     CLASS(ElemshapeData_), INTENT(IN) :: obj
     TYPE(FEVariable_), INTENT(IN) :: val
     TYPE(FEVariable_) :: ans
   END FUNCTION Interpolation1
-END INTERFACE
-
-INTERFACE Interpolation
-  MODULE PROCEDURE Interpolation1
 END INTERFACE Interpolation
 
+!----------------------------------------------------------------------------
+!
+!----------------------------------------------------------------------------
+
 END MODULE ElemshapeData_InterpolMethods
diff --git a/src/modules/FEVariable/src/FEVariable_ConstructorMethod.F90 b/src/modules/FEVariable/src/FEVariable_ConstructorMethod.F90
index 4aac79115..0bf8acba4 100644
--- a/src/modules/FEVariable/src/FEVariable_ConstructorMethod.F90
+++ b/src/modules/FEVariable/src/FEVariable_ConstructorMethod.F90
@@ -15,11 +15,9 @@
 ! along with this program.  If not, see <https: //www.gnu.org/licenses/>
 
 MODULE FEVariable_ConstructorMethod
-USE BaseType, ONLY: FEVariable_, &
-                    TypeFEVariableOpt
-
+USE BaseType, ONLY: FEVariable_
+USE BaseType, ONLY: TypeFEVariableOpt
 USE GlobalData, ONLY: I4B, DFP, LGT
-
 IMPLICIT NONE
 
 PRIVATE
diff --git a/src/modules/FEVariable/src/FEVariable_InterpolationMethod.F90 b/src/modules/FEVariable/src/FEVariable_InterpolationMethod.F90
index 1d06938b2..db53032d0 100644
--- a/src/modules/FEVariable/src/FEVariable_InterpolationMethod.F90
+++ b/src/modules/FEVariable/src/FEVariable_InterpolationMethod.F90
@@ -15,11 +15,9 @@
 ! along with this program.  If not, see <https: //www.gnu.org/licenses/>
 
 MODULE FEVariable_InterpolationMethod
-USE BaseType, ONLY: FEVariable_, &
-                    TypeFEVariableOpt
-
+USE BaseType, ONLY: FEVariable_
+USE BaseType, ONLY: TypeFEVariableOpt
 USE GlobalData, ONLY: I4B, DFP, LGT
-
 IMPLICIT NONE
 PRIVATE
 
diff --git a/src/modules/FEVariable/src/FEVariable_QuadratureVariableMethod.F90 b/src/modules/FEVariable/src/FEVariable_QuadratureVariableMethod.F90
index 4e4d3d444..a40c54f13 100644
--- a/src/modules/FEVariable/src/FEVariable_QuadratureVariableMethod.F90
+++ b/src/modules/FEVariable/src/FEVariable_QuadratureVariableMethod.F90
@@ -200,6 +200,25 @@ MODULE PURE FUNCTION Quadrature_Vector_Constant(val, rank, vartype) &
   END FUNCTION Quadrature_Vector_Constant
 END INTERFACE QuadratureVariable
 
+!----------------------------------------------------------------------------
+!                                           NodalVariable@ConstructorMethods
+!----------------------------------------------------------------------------
+
+!> author: Vikas Sharma, Ph. D.
+! date: 2021-12-10
+! update: 2021-12-10
+! summary: Create quadrature variable, which is Vector, Constant
+
+INTERFACE QuadratureVariable
+  MODULE PURE FUNCTION Quadrature_Vector_Constant2(tsize, rank, vartype) &
+    RESULT(obj)
+    TYPE(FEVariable_) :: obj
+    INTEGER(I4B), INTENT(IN) :: tsize
+    TYPE(FEVariableVector_), INTENT(IN) :: rank
+    TYPE(FEVariableConstant_), INTENT(IN) :: vartype
+  END FUNCTION Quadrature_Vector_Constant2
+END INTERFACE QuadratureVariable
+
 !----------------------------------------------------------------------------
 !                                      QuadratureVariable@ConstructorMethods
 !----------------------------------------------------------------------------
diff --git a/src/modules/GlobalData/src/GlobalData.F90 b/src/modules/GlobalData/src/GlobalData.F90
index 59cbcd4fa..17457f6d6 100755
--- a/src/modules/GlobalData/src/GlobalData.F90
+++ b/src/modules/GlobalData/src/GlobalData.F90
@@ -339,6 +339,7 @@ MODULE GlobalData
 !!
 !! Type of quadrature points
 !!
+INTEGER(I4B), PARAMETER :: CenterQP = 0
 INTEGER(I4B), PARAMETER :: EquidistanceQP = Equidistance
 INTEGER(I4B), PARAMETER :: GaussQP = Gauss
 INTEGER(I4B), PARAMETER :: GaussLegendreQP = GaussLegendre
diff --git a/src/modules/IterationData/CMakeLists.txt b/src/modules/IterationData/CMakeLists.txt
deleted file mode 100644
index e58ecfa3b..000000000
--- a/src/modules/IterationData/CMakeLists.txt
+++ /dev/null
@@ -1,13 +0,0 @@
-# This file is a part of easifem-base
-# (c) 2021, Vikas Sharma, Ph.D.
-# All right reserved
-#
-# log
-# 16/02/2021 this file was created
-#-----------------------------------------------------------------------
-
-SET(src_path "${CMAKE_CURRENT_LIST_DIR}/src/")
-TARGET_SOURCES(
-  ${PROJECT_NAME} PRIVATE
-  ${src_path}/IterationData_Method.F90
-)
\ No newline at end of file
diff --git a/src/modules/IterationData/src/IterationData_Method.F90 b/src/modules/IterationData/src/IterationData_Method.F90
deleted file mode 100644
index 34eda0561..000000000
--- a/src/modules/IterationData/src/IterationData_Method.F90
+++ /dev/null
@@ -1,108 +0,0 @@
-! This program is a part of EASIFEM library
-! Copyright (C) 2020-2021  Vikas Sharma, Ph.D
-!
-! 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 3 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.  If not, see <https: //www.gnu.org/licenses/>
-!
-
-MODULE IterationData_Method
-USE GlobalData
-USE BaseType
-IMPLICIT NONE
-PRIVATE
-
-!----------------------------------------------------------------------------
-!                                                      Initiate@Constructor
-!----------------------------------------------------------------------------
-
-INTERFACE
-MODULE PURE SUBROUTINE iterdata_Initiate( obj, maxIter, iterationNumber, &
-  & residualError0, residualError, residualTolerance, solutionError0, &
-  & solutionError, solutionTolerance, convergenceType, &
-  & convergenceIn, normType, converged, timeAtStart, timeAtEnd )
-  TYPE( IterationData_ ), INTENT( INOUT ) :: obj
-  INTEGER( I4B ), INTENT( IN ), OPTIONAL :: maxIter
-  INTEGER( I4B ), INTENT( IN ), OPTIONAL :: iterationNumber
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: residualError0
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: residualError
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: residualTolerance
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: solutionError0
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: solutionError
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: solutionTolerance
-  INTEGER( I4B ), OPTIONAL, INTENT( IN ) :: convergenceType
-  INTEGER( I4B ), OPTIONAL, INTENT( IN ) :: convergenceIn
-  INTEGER( I4B ), OPTIONAL, INTENT( IN ) :: normType
-  LOGICAL( LGT ), OPTIONAL, INTENT( IN ) :: converged
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: timeAtStart
-  REAL( DFP ), OPTIONAL, INTENT( IN ) :: timeAtEnd
-END SUBROUTINE iterdata_Initiate
-END INTERFACE
-
-INTERFACE Initiate
-  MODULE PROCEDURE iterdata_Initiate
-END INTERFACE Initiate
-
-PUBLIC :: Initiate
-
-!----------------------------------------------------------------------------
-!                                                              Deallocate
-!----------------------------------------------------------------------------
-
-INTERFACE
-MODULE SUBROUTINE iterdata_Deallocate( obj )
-  TYPE( IterationData_ ), INTENT( INOUT ) :: obj
-END SUBROUTINE iterdata_Deallocate
-END INTERFACE
-
-INTERFACE Deallocate
-  MODULE PROCEDURE iterdata_Deallocate
-END INTERFACE Deallocate
-
-PUBLIC :: Deallocate
-
-!----------------------------------------------------------------------------
-!                                                    isConverged@Constructor
-!----------------------------------------------------------------------------
-
-INTERFACE
-MODULE PURE FUNCTION iterdata_isConverged( obj ) RESULT( Ans )
-  TYPE( IterationData_ ), INTENT( IN ) :: obj
-  LOGICAL( LGT ) :: Ans
-END FUNCTION iterdata_isConverged
-END INTERFACE
-
-INTERFACE isConverged
-  MODULE PROCEDURE iterdata_isConverged
-END INTERFACE isConverged
-
-PUBLIC :: isConverged
-
-!----------------------------------------------------------------------------
-!                                                                Display@IO
-!----------------------------------------------------------------------------
-
-INTERFACE
-MODULE SUBROUTINE iterdata_Display( obj, msg, UnitNo )
-  TYPE( IterationData_ ), INTENT( IN ) :: obj
-  CHARACTER( LEN = * ), INTENT( IN ) :: msg
-  INTEGER( I4B ), OPTIONAL, INTENT( IN ) :: UnitNo
-END SUBROUTINE iterdata_Display
-END INTERFACE
-
-INTERFACE Display
-  MODULE PROCEDURE iterdata_Display
-END INTERFACE Display
-
-PUBLIC :: Display
-
-END MODULE IterationData_Method
\ No newline at end of file
diff --git a/src/modules/Lapack/src/GE_CompRoutineMethods.F90 b/src/modules/Lapack/src/GE_CompRoutineMethods.F90
index 870b57bb6..b718d08f0 100644
--- a/src/modules/Lapack/src/GE_CompRoutineMethods.F90
+++ b/src/modules/Lapack/src/GE_CompRoutineMethods.F90
@@ -27,33 +27,30 @@ MODULE GE_CompRoutineMethods
 !                                                                ConditionNo
 !----------------------------------------------------------------------------
 
-INTERFACE
+INTERFACE ConditionNo
   MODULE FUNCTION ge_ConditionNo_1(A, NORM) RESULT(ans)
     REAL(DFP), INTENT(IN) :: A(:, :)
     !! General matrix
     CHARACTER(1), INTENT(IN) :: NORM
     !! "1", "0"
     REAL(DFP) :: ans
-    !!
+    !! condition number
   END FUNCTION ge_ConditionNo_1
-END INTERFACE
-
-INTERFACE ConditionNo
-  MODULE PROCEDURE ge_ConditionNo_1
 END INTERFACE ConditionNo
 
 !----------------------------------------------------------------------------
-!                                                                    GetInvMat
+!                                                                  GetInvMat
 !----------------------------------------------------------------------------
 
 !> author: Vikas Sharma, Ph. D.
 ! date: 18 Oct 2022
-! summary:         Inverse of matrix
+! summary: Inverse of matrix
 !
-!# Introduction
+!# GetInvMat
 !
 ! This routine calls `DGETRI` routine from Lapack.
-! A copy of matrix A is made into invA, then LU decomposition is performed and
+! A copy of matrix A is made into invA,
+! then LU decomposition is performed and
 ! `DGETRI` is called from lapack
 
 INTERFACE
@@ -70,14 +67,14 @@ END SUBROUTINE ge_GetInvMat1
 END INTERFACE GetInvMat
 
 !----------------------------------------------------------------------------
-!                                                                    GetInvMat
+!                                                                  GetInvMat
 !----------------------------------------------------------------------------
 
 !> author: Vikas Sharma, Ph. D.
 ! date: 18 Oct 2022
-! summary:         Inverse of matrix
+! summary: Inverse of matrix
 !
-!# Introduction
+!# GetInvMat
 !
 !- This routine calls `DGETRI` routine from Lapack.
 !- A and IPIV are obtained from LU decomposition
@@ -85,7 +82,7 @@ END SUBROUTINE ge_GetInvMat1
 !- A copy of matrix A is made into invA, then
 ! `DGETRI` is called from lapack
 
-INTERFACE
+INTERFACE GetInvMat
   MODULE SUBROUTINE ge_GetInvMat2(A, IPIV, invA)
     REAL(DFP), INTENT(IN) :: A(:, :)
     !! General matrix
@@ -94,21 +91,17 @@ MODULE SUBROUTINE ge_GetInvMat2(A, IPIV, invA)
     REAL(DFP), INTENT(INOUT) :: invA(:, :)
     !!
   END SUBROUTINE ge_GetInvMat2
-END INTERFACE
-
-INTERFACE GetInvMat
-  MODULE PROCEDURE ge_GetInvMat2
 END INTERFACE GetInvMat
 
 !----------------------------------------------------------------------------
-!                                                                    GetInvMat
+!                                                                  GetInvMat
 !----------------------------------------------------------------------------
 
 !> author: Vikas Sharma, Ph. D.
 ! date: 18 Oct 2022
-! summary:         Inverse of matrix
+! summary: Inverse of matrix
 !
-!# Introduction
+!# GetInvMat
 !
 !- This routine calls `DGETRI` routine from Lapack.
 !- A and IPIV are obtained from LU decomposition
@@ -116,7 +109,7 @@ END SUBROUTINE ge_GetInvMat2
 !- At output invese of A is stored inside A
 !- No copy is made.
 
-INTERFACE
+INTERFACE GetInvMat
   MODULE SUBROUTINE ge_GetInvMat3(A, IPIV)
     REAL(DFP), INTENT(INOUT) :: A(:, :)
     !! LU Decompose at input
@@ -124,10 +117,6 @@ MODULE SUBROUTINE ge_GetInvMat3(A, IPIV)
     INTEGER(I4B), INTENT(IN) :: IPIV(:)
     !!
   END SUBROUTINE ge_GetInvMat3
-END INTERFACE
-
-INTERFACE GetInvMat
-  MODULE PROCEDURE ge_GetInvMat3
 END INTERFACE GetInvMat
 
 !----------------------------------------------------------------------------
@@ -136,7 +125,7 @@ END SUBROUTINE ge_GetInvMat3
 
 !> author: Vikas Sharma, Ph. D.
 ! date: 18 Oct 2022
-! summary:         Inverse of matrix
+! summary: Inverse of matrix
 !
 !# Introduction
 !
@@ -145,14 +134,10 @@ END SUBROUTINE ge_GetInvMat3
 !- Then `DGETRI` is called from lapack
 !- At output A contains the inverse.
 
-INTERFACE
+INTERFACE GetInvMat
   MODULE SUBROUTINE ge_GetInvMat4(A)
     REAL(DFP), INTENT(INOUT) :: A(:, :)
   END SUBROUTINE ge_GetInvMat4
-END INTERFACE
-
-INTERFACE GetInvMat
-  MODULE PROCEDURE ge_GetInvMat4
 END INTERFACE GetInvMat
 
 END MODULE GE_CompRoutineMethods
diff --git a/src/submodules/CMakeLists.txt b/src/submodules/CMakeLists.txt
index 264ddd354..6ce89de26 100644
--- a/src/submodules/CMakeLists.txt
+++ b/src/submodules/CMakeLists.txt
@@ -72,9 +72,6 @@ include(${CMAKE_CURRENT_LIST_DIR}/IntVector/CMakeLists.txt)
 # IndexValue
 include(${CMAKE_CURRENT_LIST_DIR}/IndexValue/CMakeLists.txt)
 
-# IterationData
-include(${CMAKE_CURRENT_LIST_DIR}/IterationData/CMakeLists.txt)
-
 # KeyValue
 include(${CMAKE_CURRENT_LIST_DIR}/KeyValue/CMakeLists.txt)
 
@@ -155,3 +152,6 @@ include(${CMAKE_CURRENT_LIST_DIR}/CSRSparsity/CMakeLists.txt)
 
 # CSRMatrix
 include(${CMAKE_CURRENT_LIST_DIR}/CSRMatrix/CMakeLists.txt)
+
+# ConvergenceOpt
+include(${CMAKE_CURRENT_LIST_DIR}/ConvergenceOpt/CMakeLists.txt)
diff --git a/src/submodules/ConvergenceOpt/CMakeLists.txt b/src/submodules/ConvergenceOpt/CMakeLists.txt
new file mode 100644
index 000000000..c3e8980c5
--- /dev/null
+++ b/src/submodules/ConvergenceOpt/CMakeLists.txt
@@ -0,0 +1,24 @@
+# This program is a part of EASIFEM library
+# Expandable And Scalable Infrastructure for Finite Element Methods
+# htttps://www.easifem.com
+# Vikas Sharma, Ph.D., vickysharma0812@gmail.com
+#
+# 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 3 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.  If not, see <https: //www.gnu.org/licenses/>
+#
+
+set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/")
+target_sources(
+    ${PROJECT_NAME}
+    PRIVATE ${src_path}/ConvergenceOptUtility@Methods.F90
+)
diff --git a/src/submodules/ConvergenceOpt/src/ConvergenceOptUtility@Methods.F90 b/src/submodules/ConvergenceOpt/src/ConvergenceOptUtility@Methods.F90
new file mode 100644
index 000000000..80f09940a
--- /dev/null
+++ b/src/submodules/ConvergenceOpt/src/ConvergenceOptUtility@Methods.F90
@@ -0,0 +1,284 @@
+! This program is a part of EASIFEM library
+! Expandable And Scalable Infrastructure for Finite Element Methods
+! htttps://www.easifem.com
+!
+! 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 3 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.  If not, see <https: //www.gnu.org/licenses/>
+!
+
+SUBMODULE(ConvergenceOptUtility) Methods
+USE Display_Method, ONLY: ToString
+USE StringUtility, ONLY: UpperCase
+USE BaseType, ONLY: TypeConvergenceOpt
+USE BaseType, ONLY: math => TypeMathOpt
+IMPLICIT NONE
+
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: modName = "ConvergenceOptUtility@Methods.F90"
+#endif
+
+CONTAINS
+
+!----------------------------------------------------------------------------
+!                                                      ConvergenceType_ToInt
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE ConvergenceType_ToInt
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: myName = "ConvergenceType_ToInt()"
+#endif
+CHARACTER(:), ALLOCATABLE :: acase
+
+acase = UpperCase(TRIM(name))
+
+SELECT CASE (acase)
+CASE ("RELATIVE")
+  ans = TypeConvergenceOpt%relative
+CASE ("ABSOLUTE")
+  ans = TypeConvergenceOpt%absolute
+CASE ("BOTH")
+  ans = TypeConvergenceOpt%both
+CASE DEFAULT
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
+                    "no case found for name="//acase)
+#endif
+END SELECT
+
+acase = ""
+END PROCEDURE ConvergenceType_ToInt
+
+!----------------------------------------------------------------------------
+!                                                      ConvergenceType_ToChar
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE ConvergenceType_ToChar
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: myName = "ConvergenceType_ToChar()"
+#endif
+
+! internal variables
+LOGICAL(LGT) :: isUpper0
+
+isUpper0 = math%no
+IF (PRESENT(isUpper)) isUpper0 = isUpper
+
+SELECT CASE (name)
+
+CASE (TypeConvergenceOpt%relative)
+  IF (isUpper0) THEN
+    ans = "RELATIVE"
+  ELSE
+    ans = "Relative"
+  END IF
+
+CASE (TypeConvergenceOpt%absolute)
+  IF (isUpper0) THEN
+    ans = "ABSOLUTE"
+  ELSE
+    ans = "Absolute"
+  END IF
+
+CASE (TypeConvergenceOpt%both)
+  IF (isUpper0) THEN
+    ans = "BOTH"
+  ELSE
+    ans = "Both"
+  END IF
+
+CASE DEFAULT
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
+                    "no case found for name="//ToString(name))
+#endif
+END SELECT
+
+END PROCEDURE ConvergenceType_ToChar
+
+!----------------------------------------------------------------------------
+!                                                   ConvergenceType_ToString
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE ConvergenceType_ToString
+ans = ConvergenceType_ToChar(name=name, isUpper=isUpper)
+END PROCEDURE ConvergenceType_ToString
+
+!----------------------------------------------------------------------------
+!                                                       ConvergenceIn_ToInt
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE ConvergenceIn_ToInt
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: myName = "ConvergenceIn_ToInt()"
+#endif
+CHARACTER(:), ALLOCATABLE :: acase
+
+acase = UpperCase(TRIM(name))
+
+SELECT CASE (acase)
+CASE ("RESIDUAL")
+  ans = TypeConvergenceOpt%relative
+CASE ("SOLUTION")
+  ans = TypeConvergenceOpt%absolute
+CASE ("BOTH")
+  ans = TypeConvergenceOpt%both
+CASE DEFAULT
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
+                    "no case found for name="//acase)
+#endif
+END SELECT
+
+acase = ""
+END PROCEDURE ConvergenceIn_ToInt
+
+!----------------------------------------------------------------------------
+!                                                        ConvergenceIn_ToChar
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE ConvergenceIn_ToChar
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: myName = "ConvergenceIn_ToChar()"
+#endif
+
+! internal variables
+LOGICAL(LGT) :: isUpper0
+
+isUpper0 = math%no
+IF (PRESENT(isUpper)) isUpper0 = isUpper
+
+SELECT CASE (name)
+
+CASE (TypeConvergenceOpt%res)
+  IF (isUpper0) THEN
+    ans = "RESIDUAL"
+  ELSE
+    ans = "Residual"
+  END IF
+
+CASE (TypeConvergenceOpt%sol)
+  IF (isUpper0) THEN
+    ans = "SOLUTION"
+  ELSE
+    ans = "Solution"
+  END IF
+
+CASE (TypeConvergenceOpt%both)
+  IF (isUpper0) THEN
+    ans = "BOTH"
+  ELSE
+    ans = "Both"
+  END IF
+
+CASE DEFAULT
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
+                    "no case found for name="//ToString(name))
+#endif
+END SELECT
+
+END PROCEDURE ConvergenceIn_ToChar
+
+!----------------------------------------------------------------------------
+!                                                     ConvergenceIn_ToString
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE ConvergenceIn_ToString
+ans = ConvergenceIn_ToChar(name=name, isUpper=isUpper)
+END PROCEDURE ConvergenceIn_ToString
+
+!----------------------------------------------------------------------------
+!                                                             NormType_ToInt
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE NormType_ToInt
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: myName = "NormType_ToInt()"
+#endif
+CHARACTER(:), ALLOCATABLE :: acase
+
+acase = UpperCase(TRIM(name))
+
+SELECT CASE (acase)
+CASE ("L1")
+  ans = TypeConvergenceOpt%normL1
+CASE ("L2")
+  ans = TypeConvergenceOpt%normL2
+CASE ("INFINITY")
+  ans = TypeConvergenceOpt%normInfinity
+
+CASE DEFAULT
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
+                    "no case found for name="//acase)
+#endif
+END SELECT
+
+acase = ""
+END PROCEDURE NormType_ToInt
+
+!----------------------------------------------------------------------------
+!                                                             NormType_ToChar
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE NormType_ToChar
+#ifdef DEBUG_VER
+CHARACTER(*), PARAMETER :: myName = "NormType_ToChar()"
+#endif
+
+! internal variables
+LOGICAL(LGT) :: isUpper0
+
+isUpper0 = math%no
+IF (PRESENT(isUpper)) isUpper0 = isUpper
+
+SELECT CASE (name)
+
+CASE (TypeConvergenceOpt%normL1)
+  ans = "L1"
+
+CASE (TypeConvergenceOpt%normL2)
+  ans = "L2"
+
+CASE (TypeConvergenceOpt%normInfinity)
+  IF (isUpper0) THEN
+    ans = "INFINITY"
+  ELSE
+    ans = "Infinity"
+  END IF
+
+CASE DEFAULT
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
+                    "no case found for name="//ToString(name))
+#endif
+END SELECT
+
+END PROCEDURE NormType_ToChar
+
+!----------------------------------------------------------------------------
+!                                                          NormType_ToString
+!----------------------------------------------------------------------------
+
+MODULE PROCEDURE NormType_ToString
+ans = NormType_ToChar(name=name, isUpper=isUpper)
+END PROCEDURE NormType_ToString
+
+!----------------------------------------------------------------------------
+!                                                              Include error
+!----------------------------------------------------------------------------
+
+#include "../../include/errors.F90"
+
+END SUBMODULE Methods
+
diff --git a/src/submodules/ElemshapeData/src/ElemshapeData_InterpolMethods@Methods.F90 b/src/submodules/ElemshapeData/src/ElemshapeData_InterpolMethods@Methods.F90
index 9f10658b5..607e2b4c0 100644
--- a/src/submodules/ElemshapeData/src/ElemshapeData_InterpolMethods@Methods.F90
+++ b/src/submodules/ElemshapeData/src/ElemshapeData_InterpolMethods@Methods.F90
@@ -17,12 +17,12 @@
 
 SUBMODULE(ElemshapeData_InterpolMethods) Methods
 USE BaseType, ONLY: TypeFEVariableOpt
-USE FEVariable_Method, ONLY: FEVariableGetInterpolation_ => GetInterpolation_,&
-                             FEVariableInitiate => Initiate, &
-                             FEVariableGetRank => GetRank, &
-                             FEVariableGetTotalShape => GetTotalShape, &
-                             FEVariableSize => Size
-
+USE BaseType, ONLY: math => TypeMathOpt
+USE FEVariable_Method, ONLY: FEVariableGetInterpolation_ => GetInterpolation_
+USE FEVariable_Method, ONLY: FEVariableInitiate => Initiate
+USE FEVariable_Method, ONLY: FEVariableGetRank => GetRank
+USE FEVariable_Method, ONLY: FEVariableGetTotalShape => GetTotalShape
+USE FEVariable_Method, ONLY: FEVariableSize => Size
 IMPLICIT NONE
 CONTAINS
 
@@ -59,6 +59,8 @@
     s(3) = obj%nips
     mylen = s(1) * s(2) * s(3)
 
+  CASE DEFAULT
+
   END SELECT
 
   CALL FEVariableInitiate(obj=ans, &
@@ -77,11 +79,9 @@
 !----------------------------------------------------------------------------
 
 MODULE PROCEDURE GetInterpolation_1
-REAL(DFP), PARAMETER :: one = 1.0_DFP
-LOGICAL, PARAMETER :: no = .FALSE.
-
 CALL FEVariableGetInterpolation_(obj=val, ans=ans, N=obj%N, nns=obj%nns, &
-                                 nips=obj%nips, scale=one, addContribution=no)
+                                 nips=obj%nips, scale=math%one, &
+                                 addContribution=math%no)
 END PROCEDURE GetInterpolation_1
 
 !----------------------------------------------------------------------------
@@ -144,6 +144,8 @@
     s(4) = nipt
     mylen = s(1) * s(2) * s(3) * s(4)
 
+  CASE DEFAULT
+
   END SELECT
 
   CALL FEVariableInitiate(obj=ans, &
@@ -158,15 +160,12 @@
 END PROCEDURE GetInterpolation2
 
 !----------------------------------------------------------------------------
-!                                                         GetInterpolation_
+!                                                           GetInterpolation_
 !----------------------------------------------------------------------------
 
 MODULE PROCEDURE GetInterpolation_2
-REAL(DFP), PARAMETER :: one = 1.0_DFP
-LOGICAL, PARAMETER :: no = .FALSE.
-
-CALL GetInterpolation_(obj=obj, ans=ans, val=val, scale=one, &
-                       addContribution=no)
+CALL GetInterpolation_(obj=obj, ans=ans, val=val, scale=math%one, &
+                       addContribution=math%no)
 END PROCEDURE GetInterpolation_2
 
 !----------------------------------------------------------------------------
@@ -188,11 +187,14 @@
 END PROCEDURE GetInterpolation_2a
 
 !----------------------------------------------------------------------------
-!                                                      interpolationOfVector
+!                                                              Interpolation
 !----------------------------------------------------------------------------
 
 MODULE PROCEDURE Interpolation1
 CALL GetInterpolation(obj=obj, val=val, ans=ans)
 END PROCEDURE Interpolation1
 
+!----------------------------------------------------------------------------
+!
+!----------------------------------------------------------------------------
 END SUBMODULE Methods
diff --git a/src/submodules/FEVariable/src/FEVariable_QuadratureVariableMethod@Methods.F90 b/src/submodules/FEVariable/src/FEVariable_QuadratureVariableMethod@Methods.F90
index d8acee37d..9a74edb6c 100644
--- a/src/submodules/FEVariable/src/FEVariable_QuadratureVariableMethod@Methods.F90
+++ b/src/submodules/FEVariable/src/FEVariable_QuadratureVariableMethod@Methods.F90
@@ -28,7 +28,6 @@
 
 MODULE PROCEDURE Quadrature_Scalar_Constant
 INTEGER(I4B) :: s(1)
-
 s(1) = 1
 CALL FEVariableInitiate(obj=obj, s=s, &
                         defineon=TypeFEVariableOpt%quadrature, &
@@ -170,6 +169,21 @@
 !                                                          QuadratureVariable
 !----------------------------------------------------------------------------
 
+MODULE PROCEDURE Quadrature_Vector_Constant2
+INTEGER(I4B) :: s(1)
+
+s(1) = tsize
+
+CALL FEVariableInitiate( &
+  obj=obj, s=s, defineon=TypeFEVariableOpt%Quadrature, &
+  vartype=TypeFEVariableOpt%constant, &
+  rank=TypeFEVariableOpt%vector, len=tsize)
+END PROCEDURE Quadrature_Vector_Constant2
+
+!----------------------------------------------------------------------------
+!                                                          QuadratureVariable
+!----------------------------------------------------------------------------
+
 MODULE PROCEDURE Quadrature_Vector_Space
 INTEGER(I4B) :: s(2), tsize, ii, jj, cnt
 
diff --git a/src/submodules/IterationData/src/IterationData_Method@IOMethods.F90 b/src/submodules/IterationData/src/IterationData_Method@IOMethods.F90
deleted file mode 100644
index d316aa439..000000000
--- a/src/submodules/IterationData/src/IterationData_Method@IOMethods.F90
+++ /dev/null
@@ -1,44 +0,0 @@
-! This program is a part of EASIFEM library
-! Copyright (C) 2020-2021  Vikas Sharma, Ph.D
-!
-! 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 3 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.  If not, see <https: //www.gnu.org/licenses/>
-!
-SUBMODULE(IterationData_Method) IOMethods
-USE BaseMethod
-IMPLICIT NONE
-CONTAINS
-
-!----------------------------------------------------------------------------
-!                                                                   Display
-!----------------------------------------------------------------------------
-
-MODULE PROCEDURE iterdata_Display
-CALL Display(TRIM(msg), unitno)
-CALL Display(obj%maxIter, 'maxIter: ', unitno)
-CALL Display(obj%iterationNumber, 'iterationNumber: ', unitno)
-CALL Display(obj%residualError0, 'residualError0: ', unitno)
-CALL Display(obj%residualError, 'residualError: ', unitno)
-CALL Display(obj%residualTolerance, 'residualTolerance: ', unitno)
-CALL Display(obj%solutionError0, 'solutionError0: ', unitno)
-CALL Display(obj%solutionError, 'solutionError: ', unitno)
-CALL Display(obj%solutionTolerance, 'solutionTolerance: ', unitno)
-CALL Display(obj%convergenceType, 'convergenceType: ', unitno)
-CALL Display(obj%convergenceIn, 'convergenceIn: ', unitno)
-CALL Display(obj%normType, 'normType: ', unitno)
-CALL Display(obj%converged, 'converged: ', unitno)
-CALL Display(obj%timeAtStart, 'timeAtStart: ', unitno)
-CALL Display(obj%timeAtEnd, 'timeAtEnd: ', unitno)
-END PROCEDURE iterdata_Display
-
-END SUBMODULE IOMethods
diff --git a/src/submodules/Triangle/src/TriangleInterpolationUtility@Methods.F90 b/src/submodules/Triangle/src/TriangleInterpolationUtility@Methods.F90
index e5119a32b..f2971e160 100644
--- a/src/submodules/Triangle/src/TriangleInterpolationUtility@Methods.F90
+++ b/src/submodules/Triangle/src/TriangleInterpolationUtility@Methods.F90
@@ -16,17 +16,18 @@
 
 SUBMODULE(TriangleInterpolationUtility) Methods
 USE BaseType, ONLY: ipopt => TypeInterpolationOpt
+USE BaseType, ONLY: math => TypeMathOpt
 USE StringUtility, ONLY: UpperCase
 USE MappingUtility, ONLY: FromUnitTriangle2Triangle_
 USE RecursiveNodesUtility, ONLY: RecursiveNode2D_
 USE Display_Method, ONLY: ToString
 USE IntegerUtility, ONLY: NumberOfTuples => SIZE
-USE LineInterpolationUtility, ONLY: EquidistanceInPoint_Line_, &
-                                    InterpolationPoint_Line_
-
+USE LineInterpolationUtility, ONLY: EquidistanceInPoint_Line_
+USE LineInterpolationUtility, ONLY: InterpolationPoint_Line_
 IMPLICIT NONE
 
-CHARACTER(*), PARAMETER :: modName = "TriangleInterpolationUtility%Methods"
+CHARACTER(*), PARAMETER :: modName = &
+                           "TriangleInterpolationUtility@Methods.F90"
 
 CONTAINS
 
@@ -79,20 +80,24 @@
   CASE ("O")
     ans = "BIUNIT"
 
-#ifdef DEBUG_VER
   CASE DEFAULT
-    CALL AssertError1(.FALSE., myName, modName, __LINE__, &
-                  "No case found for given baseInterpol="//TRIM(baseInterpol))
+
+#ifdef DEBUG_VER
+    CALL AssertError1(math%no, myName, modName, &
+                      __LINE__, &
+                      "No case found for given baseInterpol="// &
+                      TRIM(baseInterpol))
 #endif
 
   END SELECT
 
-#ifdef DEBUG_VER
 CASE DEFAULT
-  CALL AssertError1(.FALSE., myName, modName, __LINE__, &
-              "No case found for given baseContinuity="//TRIM(baseContinuity))
+#ifdef DEBUG_VER
+  CALL AssertError1(math%no, myName, modName, &
+                    __LINE__, &
+                    "No case found for given baseContinuity="// &
+                    TRIM(baseContinuity))
 #endif
-
 END SELECT
 
 END PROCEDURE RefElemDomain_Triangle
@@ -102,29 +107,9 @@
 !----------------------------------------------------------------------------
 
 MODULE PROCEDURE FacetConnectivity_Triangle
-! CHARACTER(1) :: bi
-! LOGICAL(LGT) :: isok
-
 ans(1:2, 1) = [1, 2]
 ans(1:2, 2) = [2, 3]
 ans(1:2, 3) = [3, 1]
-
-! isok = PRESENT(baseInterpol)
-! bi = "L"
-! IF (isok) bi = UpperCase(baseInterpol(1:1))
-!
-! SELECT CASE (bi)
-! CASE ("H", "O")
-!   ans(1:2, 1) = [1, 2]
-!   ans(1:2, 2) = [1, 3]
-!   ans(1:2, 3) = [2, 3]
-!
-! CASE DEFAULT
-!   ans(1:2, 1) = [1, 2]
-!   ans(1:2, 2) = [2, 3]
-!   ans(1:2, 3) = [3, 1]
-!
-! END SELECT
 END PROCEDURE FacetConnectivity_Triangle
 
 !----------------------------------------------------------------------------
@@ -155,7 +140,10 @@
 REAL(DFP) :: x(3, 3), xin(3, 3), e1(3), e2(3), lam, avar, mu
 LOGICAL(LGT) :: isok
 
-x = 0.0_DFP; xin = 0.0_DFP; e1 = 0.0_DFP; e2 = 0.0_DFP
+x = 0.0_DFP
+xin = 0.0_DFP
+e1 = 0.0_DFP
+e2 = 0.0_DFP
 
 isok = PRESENT(xij)
 
@@ -164,9 +152,9 @@
   x(1:nrow, 1:3) = xij(1:nrow, 1:3)
 ELSE
   nrow = 2_I4B
-  x(1:nrow, 1) = [0.0, 0.0]
-  x(1:nrow, 2) = [1.0, 0.0]
-  x(1:nrow, 3) = [0.0, 1.0]
+  x(1:nrow, 1) = [math%zero, math%zero]
+  x(1:nrow, 2) = [math%one, math%zero]
+  x(1:nrow, 3) = [math%zero, math%one]
 END IF
 
 ncol = LagrangeDOF_Triangle(order=order)
@@ -235,8 +223,10 @@
 xin(1:nrow, 3) = x(1:nrow, 3) + lam * e1(1:nrow) + mu * e2(1:nrow)
 
 i1 = i2 + 1
-CALL EquidistancePoint_Triangle_(order=order - 3, xij=xin(1:nrow, 1:3), &
-                                 ans=ans(1:nrow, i1:), nrow=aint, ncol=bint)
+CALL EquidistancePoint_Triangle_(order=order - 3, &
+                                 xij=xin(1:nrow, 1:3), &
+                                 ans=ans(1:nrow, i1:), &
+                                 nrow=aint, ncol=bint)
 END PROCEDURE EquidistancePoint_Triangle_
 
 !----------------------------------------------------------------------------
@@ -254,7 +244,8 @@
 END IF
 
 isok = PRESENT(xij)
-nrow = 2_I4B; IF (isok) nrow = SIZE(xij, 1)
+nrow = 2_I4B
+IF (isok) nrow = SIZE(xij, 1)
 ncol = LagrangeInDOF_Triangle(order=order)
 
 CALL EquidistanceInPoint_Triangle_(order=order, ans=ans, nrow=nrow, &
@@ -270,18 +261,23 @@
 INTEGER(I4B) :: aint, bint
 LOGICAL(LGT) :: isok
 
-nrow = 0; ncol = 0
+nrow = 0
+ncol = 0
 
 isok = order .LT. 3_I4B
 IF (isok) RETURN
 
-x = 0.0_DFP; xin = 0.0_DFP; e1 = 0.0_DFP; e2 = 0.0_DFP
+x = 0.0_DFP
+xin = 0.0_DFP
+e1 = 0.0_DFP
+e2 = 0.0_DFP
 
 isok = PRESENT(xij)
 nrow = 2_I4B
-x(1:nrow, 1) = [0.0, 0.0]
-x(1:nrow, 2) = [1.0, 0.0]
-x(1:nrow, 3) = [0.0, 1.0]
+x(1:nrow, 1) = [math%zero, math%zero]
+x(1:nrow, 2) = [math%one, math%zero]
+x(1:nrow, 3) = [math%zero, math%one]
+
 IF (isok) THEN
   nrow = SIZE(xij, 1)
   x(1:nrow, 1:3) = xij(1:nrow, 1:3)
@@ -341,11 +337,12 @@
 
 isok = PRESENT(xij)
 ncol = LagrangeDOF_Triangle(order=order)
-nrow = 2; IF (isok) nrow = SIZE(xij, 1)
+nrow = 2
+IF (isok) nrow = SIZE(xij, 1)
 ALLOCATE (ans(nrow, ncol))
 CALL BlythPozrikidis_Triangle_( &
-  order=order, ipType=ipType, ans=ans, nrow=nrow, ncol=ncol, layout=layout, &
-  xij=xij, alpha=alpha, beta=beta, lambda=lambda)
+  order=order, ipType=ipType, ans=ans, nrow=nrow, ncol=ncol, &
+  layout=layout, xij=xij, alpha=alpha, beta=beta, lambda=lambda)
 END PROCEDURE BlythPozrikidis_Triangle
 
 !----------------------------------------------------------------------------
@@ -374,17 +371,18 @@
 ncol = LagrangeDOF_Triangle(order=order)
 nrow = 2
 
-isx = .FALSE.; IF (PRESENT(xij)) isx = .TRUE.
+isx = math%no
+IF (PRESENT(xij)) isx = math%yes
 IF (isx) nrow = SIZE(xij, 1)
 
-xi(1:order + 1, 1:order + 1) = 0.0_DFP
-eta(1:order + 1, 1:order + 1) = 0.0_DFP
+xi(1:order + 1, 1:order + 1) = math%zero
+eta(1:order + 1, 1:order + 1) = math%zero
 
 DO ii = 1, order + 1
   DO jj = 1, order + 2 - ii
     kk = order + 3 - ii - jj
-    xi(ii, jj) = (1.0 + 2.0 * v(ii) - v(jj) - v(kk)) / 3.0_DFP
-    eta(ii, jj) = (1.0 + 2.0 * v(jj) - v(ii) - v(kk)) / 3.0_DFP
+    xi(ii, jj) = (math%one + math%two * v(ii) - v(jj) - v(kk)) / math%three
+    eta(ii, jj) = (math%one + math%two * v(jj) - v(ii) - v(kk)) / math%three
   END DO
 END DO
 
@@ -403,14 +401,15 @@
 
   ans(1:2, 1:ncol) = temp(1:2, 1:ncol)
 
-#ifdef DEBUG_VER
 CASE DEFAULT
-  CALL AssertError1(.FALSE., myName, modName, __LINE__, &
+
+#ifdef DEBUG_VER
+  CALL AssertError1(math%yes, myName, modName, &
+                    __LINE__, &
                     "layout=VEFC is allowed, found layout is "//TRIM(layout))
 #endif
 
 END SELECT
-
 END PROCEDURE BlythPozrikidis_Triangle_
 
 !----------------------------------------------------------------------------
@@ -421,14 +420,14 @@
 INTEGER(I4B) :: nrow, ncol
 
 ncol = NumberOfTuples(n=order, d=2)
-nrow = 2; IF (PRESENT(xij)) nrow = SIZE(xij, 1)
+nrow = 2
+IF (PRESENT(xij)) nrow = SIZE(xij, 1)
 
 ALLOCATE (ans(nrow, ncol))
 
 CALL Isaac_Triangle_(order=order, ipType=ipType, ans=ans, nrow=nrow, &
                      ncol=ncol, layout=layout, xij=xij, alpha=alpha, &
                      beta=beta, lambda=lambda)
-
 END PROCEDURE Isaac_Triangle
 
 !----------------------------------------------------------------------------
@@ -445,8 +444,6 @@
 REAL(DFP) :: xi(max_order + 1, max_order + 1), &
              eta(max_order + 1, max_order + 1), &
              temp(2, 512)
-
-! REAL(DFP), ALLOCATABLE :: temp(:, :), rPoints(:, :)
 INTEGER(I4B) :: cnt, ii, jj
 INTEGER(I4B) :: nn
 
@@ -461,8 +458,8 @@
 
 !! convert from rPoints to xi and eta
 cnt = 0
-xi(1:nn, 1:nn) = 0.0_DFP
-eta(1:nn, 1:nn) = 0.0_DFP
+xi(1:nn, 1:nn) = math%zero
+eta(1:nn, 1:nn) = math%zero
 
 DO ii = 1, nn
   DO jj = 1, nn + 1 - ii
@@ -488,7 +485,7 @@
 END IF
 
 #ifdef DEBUG_VER
-CALL AssertError1(.FALSE., myName, modName, __LINE__, &
+CALL AssertError1(math%no, myName, modName, __LINE__, &
                   "Only layout=VEFC is allowed, found layout is "//layout)
 #endif
 END PROCEDURE Isaac_Triangle_
@@ -512,17 +509,20 @@
 DO ll = 0, llt
   !! v1
   cnt = cnt + 1
-  ii = 1 + ll; jj = 1 + ll
+  ii = 1 + ll
+  jj = 1 + ll
   temp(1, cnt) = xi(ii, jj)
   temp(2, cnt) = eta(ii, jj)
   !! v2
   cnt = cnt + 1
-  ii = m + 1 - 2 * ll; jj = 1 + ll
+  ii = m + 1 - 2 * ll
+  jj = 1 + ll
   temp(1, cnt) = xi(ii, jj)
   temp(2, cnt) = eta(ii, jj)
   !! v3
   cnt = cnt + 1
-  ii = 1 + ll; jj = m + 1 - 2 * ll
+  ii = 1 + ll
+  jj = m + 1 - 2 * ll
   temp(1, cnt) = xi(ii, jj)
   temp(2, cnt) = eta(ii, jj)
   !! nodes on edge 12
@@ -554,7 +554,8 @@
   !! a internal point
   cnt = cnt + 1
   ll = llt + 1
-  ii = 1 + ll; jj = 1 + ll
+  ii = 1 + ll
+  jj = 1 + ll
   temp(1, cnt) = xi(ii, jj)
   temp(2, cnt) = eta(ii, jj)
 END IF
@@ -589,13 +590,19 @@
       ipopt%GaussLegendreLobatto, ipopt%GaussChebyshevLobatto)
   ncol = NumberOfTuples(n=order, d=2)
 
+CASE (ipopt%Center)
+  ncol = 1_DFP
+
+CASE DEFAULT
+
 END SELECT
 
 ALLOCATE (ans(nrow, ncol))
 
-CALL InterpolationPoint_Triangle_(order=order, ipType=ipType, ans=ans, &
-                                 nrow=nrow, ncol=ncol, xij=xij, alpha=alpha, &
-                                  beta=beta, lambda=lambda, layout=layout)
+CALL InterpolationPoint_Triangle_( &
+  order=order, ipType=ipType, ans=ans, &
+  nrow=nrow, ncol=ncol, xij=xij, alpha=alpha, &
+  beta=beta, lambda=lambda, layout=layout)
 
 END PROCEDURE InterpolationPoint_Triangle
 
@@ -609,6 +616,19 @@
 #endif
 
 SELECT CASE (ipType)
+CASE (ipopt%Center)
+  IF (PRESENT(xij)) THEN
+    nrow = SIZE(xij, 1)
+    ncol = 1
+    ans(1:nrow, ncol) = (xij(1:nrow, 1) + xij(1:nrow, 2) &
+                         + xij(1:nrow, 3)) / math%three
+  ELSE
+    nrow = 2
+    ncol = 1
+    ans(1, 1) = math%one / math%three
+    ans(2, 1) = math%one / math%three
+  END IF
+
 CASE (ipopt%Equidistance)
   CALL EquidistancePoint_Triangle_(xij=xij, order=order, ans=ans, &
                                    nrow=nrow, ncol=ncol)
@@ -640,7 +660,7 @@
 CASE (ipopt%Feket, ipopt%Hesthaven, ipopt%ChenBabuska)
 
 #ifdef DEBUG_VER
-  CALL AssertError1(.FALSE., myName, modName, __LINE__, &
+  CALL AssertError1(math%no, myName, modName, __LINE__, &
                     "Feket, Hesthaven, ChenBabuska nodes not available")
 #endif