diff --git a/src/modules/STVectorField/src/STVectorField_Class.F90 b/src/modules/STVectorField/src/STVectorField_Class.F90 index 339467bb..8a7593f3 100644 --- a/src/modules/STVectorField/src/STVectorField_Class.F90 +++ b/src/modules/STVectorField/src/STVectorField_Class.F90 @@ -19,16 +19,17 @@ ! summary: STVector field data type is defined MODULE STVectorField_Class -USE GlobalData, ONLY: DFP, I4B, LGT, DOF_FMT, NODES_FMT, NodesToDOF -USE BaseType, ONLY: FEVariable_ USE AbstractField_Class, ONLY: AbstractField_ USE AbstractNodeField_Class, ONLY: AbstractNodeField_ +USE BaseType, ONLY: FEVariable_ +USE DirichletBC_Class, ONLY: DirichletBC_, DirichletBCPointer_ USE ExceptionHandler_Class, ONLY: e -USE HDF5File_Class, ONLY: HDF5File_ USE FEDOF_Class, ONLY: FEDOF_, FEDOFPointer_ -USE DirichletBC_Class, ONLY: DirichletBC_, DirichletBCPointer_ -USE UserFunction_Class, ONLY: UserFunction_ +USE GlobalData, ONLY: DFP, I4B, LGT, DOF_FMT, NODES_FMT, NodesToDOF +USE HDF5File_Class, ONLY: HDF5File_ USE TimeFEDOF_Class, ONLY: TimeFEDOF_, TimeFEDOFPointer_ +USE UserFunction_Class, ONLY: UserFunction_ +USE VectorField_Class, ONLY: VectorField_ IMPLICIT NONE @@ -68,7 +69,10 @@ MODULE STVectorField_Class ! CONSTRUCTOR: ! @ConstructorMethods - PROCEDURE, NON_OVERRIDABLE, PUBLIC, PASS(obj) :: Initiate2 => obj_Initiate2 + PROCEDURE, PUBLIC, PASS(obj) :: Initiate2 => obj_Initiate2 + !! Initiate by copy + + PROCEDURE, PUBLIC, PASS(obj) :: Initiate4 => obj_Initiate4 !! Initiate by copy PROCEDURE, PUBLIC, PASS(obj) :: DEALLOCATE => obj_Deallocate @@ -146,10 +150,11 @@ MODULE STVectorField_Class Set7, Set8, Set9, Set10, Set11, Set12, Set13, Set14, Set15, Set16 PROCEDURE, PUBLIC, NON_OVERRIDABLE, PASS(obj) :: SetByFunction => obj_SetByFunction - + !! Set by function PROCEDURE, PUBLIC, NON_OVERRIDABLE, PASS(obj) :: SetFromVectorField => & obj_SetFromVectorField - !! Set by function + PROCEDURE, PUBLIC, NON_OVERRIDABLE, PASS(obj) :: SetToVectorField => & + obj_SetToVectorField ! GET: ! @GetMethods @@ -204,6 +209,36 @@ MODULE STVectorField_Class GENERIC, PUBLIC :: ApplyDirichletBC => ApplyDirichletBC1, & ApplyDirichletBC2, ApplyDirichletBC3 + ! SET: + ! @BodySourceMethods + PROCEDURE, NON_OVERRIDABLE, PASS(obj) :: & + ApplyBodySource1 => obj_ApplyBodySource1 + !! Add contribution of body source to the scalar field + !! body source is given as user function + PROCEDURE, NON_OVERRIDABLE, PASS(obj) :: & + ApplyBodySource2 => obj_ApplyBodySource2 + !! Add contribution of body source to the scalar field + !! body source is given external scalar field + PROCEDURE, NON_OVERRIDABLE, PASS(obj) :: & + ApplyBodySource3 => obj_ApplyBodySource3 + !! Add contribution of body source to the scalar field + !! body source is given as user function + GENERIC, PUBLIC :: ApplyBodySource => ApplyBodySource1, & + ApplyBodySource2, ApplyBodySource3 + !! Generic method for setting body source + + ! SET: + ! @SurfaceNBCMethods + PROCEDURE, PUBLIC, NON_OVERRIDABLE, PASS(obj) :: ApplySurfaceNeumannBC => & + obj_ApplySurfaceNeumannBC + !! Apply Surface neumann boundary condition + + ! SET: + ! @PointNBCMethods + PROCEDURE, PUBLIC, NON_OVERRIDABLE, PASS(obj) :: ApplyPointNeumannBC => & + obj_ApplyPointNeumannBC + !! Apply point neumann boundary condition + END TYPE STVectorField_ !--------------------------------------------------------------------------- @@ -222,7 +257,7 @@ MODULE STVectorField_Class ! date: 2023-03-29 ! summary: Initiate2 -INTERFACE STVectorFieldInitiate +INTERFACE MODULE SUBROUTINE obj_Initiate2( & obj, obj2, copyFull, copyStructure, usePointer) CLASS(STVectorField_), INTENT(INOUT) :: obj @@ -232,6 +267,93 @@ MODULE SUBROUTINE obj_Initiate2( & LOGICAL(LGT), OPTIONAL, INTENT(IN) :: copyStructure LOGICAL(LGT), OPTIONAL, INTENT(IN) :: usePointer END SUBROUTINE obj_Initiate2 +END INTERFACE + +INTERFACE STVectorFieldInitiate + MODULE PROCEDURE obj_Initiate2 +END INTERFACE STVectorFieldInitiate + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +INTERFACE + MODULE SUBROUTINE obj_Initiate4( & + obj, name, engine, fieldType, storageFMT, comm, local_n, global_n, & + spaceCompo, isSpaceCompo, isSpaceCompoScalar, timeCompo, isTimeCompo, & + isTimeCompoScalar, tPhysicalVarNames, physicalVarNames, & + isPhysicalVarNames, isPhysicalVarNamesScalar, tNodes, isTNodes, & + isTNodesScalar, tSize, fedof, geofedof, timefedof) + CLASS(STVectorField_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: name + !! name of the field, needed + CHARACTER(*), INTENT(IN) :: engine + !! name of the engine, needed + INTEGER(I4B), OPTIONAL, INTENT(IN) :: fieldType + !! field type, default is FIELD_TYPE_NORMAL, needed + !! following options are available FIELD_TYPE_NORMAL + !! FIELD_TYPE_CONSTANT + INTEGER(I4B), OPTIONAL, INTENT(IN) :: storageFMT + !! storage format of the Vector field + !! Not required. + INTEGER(I4B), OPTIONAL, INTENT(IN) :: comm + !! communication group + !! Only needed for parallel environment + INTEGER(I4B), OPTIONAL, INTENT(IN) :: local_n + !! local size of field on each processor + !! Only needed for parallel environment + INTEGER(I4B), OPTIONAL, INTENT(IN) :: global_n + !! global size of field on distributed on processors + !! Only needed for parallel environment + INTEGER(I4B), OPTIONAL, INTENT(IN) :: spaceCompo(:) + !! space components + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isSpaceCompo + !! if true we will try to access spaceCompo, NOT REQUIRED + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isSpaceCompoScalar + !! is space component scalar, + !! in this case we only access spaceCompo(1), NOT REQUIRED + INTEGER(I4B), OPTIONAL, INTENT(IN) :: timeCompo(:) + !! Time components, needed for space-time field + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isTimeCompo + !! if true we will try to access TimeCompo + !! Not required, Not needed + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isTimeCompoScalar + !! is Time component scalar, + !! in this case we only access TimeCompo(1), Not needed + INTEGER(I4B), OPTIONAL, INTENT(IN) :: tPhysicalVarNames + !! total physical variable names + !! if it is zero, then physicalVarNames will not be written + !! evenif physicalVarNames is present, and isPhysicalVarNames + !! is true, Not required + CHARACTER(*), OPTIONAL, INTENT(IN) :: physicalVarNames(:) + !! Names of the physical variables, NOT REQUIRED + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isPhysicalVarNames + !! logical variable to check if physicalVarNames is present or not + !! if it is false then physicalVarNames will not be written + !! NOT REQUIRED + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isPhysicalVarNamesScalar + !! if true then physicalVarNames is scalar + INTEGER(I4B), OPTIONAL, INTENT(IN) :: tNodes(:) + !! total number of nodes in each physical variable + !! Not required + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isTNodes + !! if true we will try to access tNodes + !! Not required + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: isTNodesScalar + !! is tNodes scalar + !! Not required + INTEGER(I4B), OPTIONAL, INTENT(IN) :: tSize + !! total size of node field + !! not required + CLASS(FEDOF_), TARGET, INTENT(IN) :: fedof, geofedof + !! FEDOF object + CLASS(TimeFEDOF_), OPTIONAL, TARGET, INTENT(IN) :: timefedof + !! TimeFEDOF object + END SUBROUTINE obj_Initiate4 +END INTERFACE + +INTERFACE STVectorFieldInitiate + MODULE PROCEDURE obj_Initiate4 END INTERFACE STVectorFieldInitiate !---------------------------------------------------------------------------- @@ -879,6 +1001,25 @@ MODULE SUBROUTINE obj_SetFromVectorField(obj, VALUE, timeCompo, & END SUBROUTINE obj_SetFromVectorField END INTERFACE +!---------------------------------------------------------------------------- +! Set@SetMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Set to VectorField_ + +INTERFACE + MODULE SUBROUTINE obj_SetToVectorField(obj, VALUE, timeCompo, & + scale, addContribution) + CLASS(STVectorField_), INTENT(INOUT) :: obj + CLASS(VectorField_), INTENT(INOUT) :: VALUE + INTEGER(I4B), INTENT(IN) :: timeCompo + REAL(DFP), OPTIONAL, INTENT(IN) :: scale + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: addContribution + END SUBROUTINE obj_SetToVectorField +END INTERFACE + !---------------------------------------------------------------------------- ! Set@SetMethods !---------------------------------------------------------------------------- @@ -1262,6 +1403,114 @@ MODULE SUBROUTINE obj_ApplyDirichletBC3(obj, times) END SUBROUTINE obj_ApplyDirichletBC3 END INTERFACE +!---------------------------------------------------------------------------- +! ApplyBodySource@NBCMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Apply Body source to ST vector field + +INTERFACE + MODULE SUBROUTINE obj_ApplyBodySource1(obj, bodySource, scale, times) + CLASS(STVectorField_), INTENT(INOUT) :: obj + CLASS(UserFunction_), INTENT(INOUT) :: bodySource + !! Body source user function + !! It should be a vector function with + !! total arguments 4 (x, y, z, time) + REAL(DFP), INTENT(IN) :: scale + !! scale for body source + !! obj = obj + scale * bodySource integral + REAL(DFP), INTENT(IN) :: times(:) + !! time, which will be passed to the body source function + END SUBROUTINE obj_ApplyBodySource1 +END INTERFACE + +!---------------------------------------------------------------------------- +! ApplyBodySource@NBCMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Apply Body source to ST vector field + +INTERFACE + MODULE SUBROUTINE obj_ApplyBodySource2(obj, bodySource, scale) + CLASS(STVectorField_), INTENT(INOUT) :: obj + !! ST vector field + !! The test function will be corresponding to the obj + CLASS(STVectorField_), INTENT(INOUT) :: bodySource + !! Body source in terms of vector field + REAL(DFP), INTENT(IN) :: scale + !! scale for body source + !! obj = obj + scale * bodySource integral + END SUBROUTINE obj_ApplyBodySource2 +END INTERFACE + +!---------------------------------------------------------------------------- +! ApplyBodySource@NBCMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Apply Body source to ST vector field + +INTERFACE + MODULE SUBROUTINE obj_ApplyBodySource3(obj, bodySource, scale, times) + CLASS(STVectorField_), INTENT(INOUT) :: obj + CLASS(UserFunction_), INTENT(INOUT) :: bodySource + !! Body source user function + !! It should be a vector function with + !! total arguments 4 (x, y, z, time) + REAL(DFP), INTENT(IN) :: scale + !! scale for body source + !! obj = obj + scale * bodySource integral + REAL(DFP), INTENT(IN) :: times + !! time, which will be passed to the body source function + !! This time can also represent a quadrature point time + END SUBROUTINE obj_ApplyBodySource3 +END INTERFACE + +!---------------------------------------------------------------------------- +! ApplyNeumannBC@NBCMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Add Contribution of neumann boundary condition + +INTERFACE + MODULE SUBROUTINE obj_ApplySurfaceNeumannBC(obj, nbcField, scale, times) + CLASS(STVectorField_), INTENT(INOUT) :: obj + !! Vector field + CLASS(STVectorField_), INTENT(INOUT) :: nbcField + !! Vector field where we will keep the neumann boundary condition + !! extension to the entire domain + REAL(DFP), INTENT(IN) :: scale + !! Scale for neumann boundary condition + REAL(DFP), INTENT(IN) :: times(:) + !! times + END SUBROUTINE obj_ApplySurfaceNeumannBC +END INTERFACE + +!---------------------------------------------------------------------------- +! ApplyPointNeumannBC@NBCMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Add Contribution of point neumann boundary condition + +INTERFACE + MODULE SUBROUTINE obj_ApplyPointNeumannBC(obj, scale, times) + CLASS(STVectorField_), INTENT(INOUT) :: obj + !! Vector field + REAL(DFP), INTENT(IN) :: scale + !! scale for neumann boundary condition + REAL(DFP), INTENT(IN) :: times(:) + END SUBROUTINE obj_ApplyPointNeumannBC +END INTERFACE + !---------------------------------------------------------------------------- ! !---------------------------------------------------------------------------- diff --git a/src/modules/TDGAlgorithms/CMakeLists.txt b/src/modules/TDGAlgorithms/CMakeLists.txt index c9b1dc5d..f739e928 100644 --- a/src/modules/TDGAlgorithms/CMakeLists.txt +++ b/src/modules/TDGAlgorithms/CMakeLists.txt @@ -16,4 +16,9 @@ # set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/") -target_sources(${PROJECT_NAME} PRIVATE ${src_path}/TDGAlgorithm1_Class.F90) +target_sources( + ${PROJECT_NAME} + PRIVATE + ${src_path}/TDGAlgorithm1_Class.F90 + ${src_path}/TDGAlgorithm2_Class.F90 +) diff --git a/src/modules/TDGAlgorithms/src/TDGAlgorithm2_Class.F90 b/src/modules/TDGAlgorithms/src/TDGAlgorithm2_Class.F90 new file mode 100644 index 00000000..7ad229e5 --- /dev/null +++ b/src/modules/TDGAlgorithms/src/TDGAlgorithm2_Class.F90 @@ -0,0 +1,267 @@ +! 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 + +MODULE TDGAlgorithm2_Class +USE GlobalData, ONLY: I4B, DFP, LGT +USE TxtFile_Class, ONLY: TxtFile_ +USE ExceptionHandler_Class, ONLY: e +USE tomlf, ONLY: toml_table +USE BaseType, ONLY: ElemShapeData_ + +IMPLICIT NONE + +PRIVATE +PUBLIC :: TDGAlgorithm2_ +CHARACTER(*), PARAMETER :: modName = "TDGAlgorithm2_Class()" + +INTEGER(I4B), PARAMETER :: MAX_ORDER_TIME = 20 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-12 +! summary: Velocity based time discontinuous Galerkin algorithm + +TYPE :: TDGAlgorithm2_ + LOGICAL(LGT) :: isInit = .FALSE. + !! Flag to check if the object is initiated + + CHARACTER(4) :: name = "TDG2" + + INTEGER(I4B) :: nrow = 0_I4B, ncol = 0_I4B + !! Number of rows and columns in ct, mt, mtplus matrices + + REAL(DFP) :: initialGuess(MAX_ORDER_TIME + 4) = 0.0_DFP + !! coefficient for initial guess of solution + LOGICAL(LGT) :: initialGuess_zero(MAX_ORDER_TIME + 4) = .TRUE. + + REAL(DFP) :: dis(MAX_ORDER_TIME + 4) = 0.0_DFP + !! dis coefficient for displacement update + !! displacement = dis(1)*Un+dis(2)*Vn*dt +dis(3)*An*dt^2+ dis(4)*sol(1) * dt + ... + !! dis(1) coefficient of displacement at time tn + !! dis(2) coefficient of velocity at time tn + !! dis(3) coefficient of acceleration at time tn + !! dis(4:MAX_ORDER_TIME+4) coefficient of solution dof at time t1, t2, ... + LOGICAL(LGT) :: dis_zero(MAX_ORDER_TIME + 4) = .TRUE. + + REAL(DFP) :: vel(MAX_ORDER_TIME + 4) = 0.0_DFP + !! vel coefficient for velocity update + !! velocity = vel(1)*Un / dt + vel(2) * Vn + vel(3)*An *dt + vel(4) * sol(1) + ... + !! vel(1) coefficient of displacement at time tn + !! vel(2) coefficient of velocity at time tn + !! vel(3) coefficient of acceleration at time tn + !! vel(4:MAX_ORDER_TIME+4) coefficient of solution dof at time t1, t2, ... + LOGICAL(LGT) :: vel_zero(MAX_ORDER_TIME + 4) = .TRUE. + + REAL(DFP) :: acc(MAX_ORDER_TIME + 4) = 0.0_DFP + !! acc coefficient for acceleration update + !! acceleration = (1)*Un / dt^2 + acc(2) * Vn^2 + acc(3) * An + acc(4) * sol(1) / dt + ... + !! acc(1) coefficient of displacement at time tn + !! acc(2) coefficient of velocity at time tn + !! acc(3) coefficient of acceleration at time tn + !! acc(4:MAX_ORDER_TIME+4) coefficient of solution dof at time t1, t2, ... + LOGICAL(LGT) :: acc_zero(MAX_ORDER_TIME + 4) = .TRUE. + + REAL(DFP) :: mt(MAX_ORDER_TIME + 1, MAX_ORDER_TIME + 1) = 0.0_DFP + !! coefficient for mass matrix M + !! This is equivalent to temporal convective matrix + + !! jump contribution (mtplus) + + REAL(DFP) :: ct(MAX_ORDER_TIME + 1, MAX_ORDER_TIME + 1) = 0.0_DFP + !! coefficient for damping matrix C*dt + !! This is equivalent to temporal mass matrix + + REAL(DFP) :: kt(MAX_ORDER_TIME + 1, MAX_ORDER_TIME + 1) = 0.0_DFP + !! coefficient for stiffness matrix K * dt^2 + !! This is equivalent to T^t*\tilde{T} + !! \tilde{T} is displacement-velocity function derived from + !! the kinematical condition \dot{u}=v + + REAL(DFP) :: bt(MAX_ORDER_TIME + 1, 2 * MAX_ORDER_TIME + 2) = 0.0_DFP + !! sub matrix used to derive kt + + REAL(DFP) :: bt_right(MAX_ORDER_TIME + 1) = 0.0_DFP + !! bt at theta +1 + + REAL(DFP) :: wt(MAX_ORDER_TIME + 1, MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: wmt(MAX_ORDER_TIME + 1, MAX_ORDER_TIME + 1) = 0.0_DFP + !! transpose of wt*mt + !! needed to make bt + !! this matrix is made when we call GetWt + + REAL(DFP) :: at(MAX_ORDER_TIME + 1) = 0.0_DFP + !! At matrix + + REAL(DFP) :: at_right = 0.0_DFP + !! At at theta +1 + + REAL(DFP) :: tat(MAX_ORDER_TIME + 1) = 0.0_DFP + !! integral of shape function of time times at + + REAL(DFP) :: rhs_m_u1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_m_v1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_m_a1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_k_u1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_k_v1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_k_a1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_c_u1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_c_v1(MAX_ORDER_TIME + 1) = 0.0_DFP + REAL(DFP) :: rhs_c_a1(MAX_ORDER_TIME + 1) = 0.0_DFP + + LOGICAL(LGT) :: rhs_m_u1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_m_v1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_m_a1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_k_u1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_k_v1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_k_a1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_c_u1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_c_v1_zero(MAX_ORDER_TIME + 1) = .TRUE. + LOGICAL(LGT) :: rhs_c_a1_zero(MAX_ORDER_TIME + 1) = .TRUE. + +CONTAINS + PRIVATE + PROCEDURE, PUBLIC, PASS(obj) :: IsInitiated => obj_IsInitiated + PROCEDURE, PUBLIC, PASS(obj) :: Initiate => obj_Initiate + PROCEDURE, PUBLIC, PASS(obj) :: DEALLOCATE => obj_Deallocate + PROCEDURE, PUBLIC, PASS(obj) :: Display => obj_Display + PROCEDURE, PUBLIC, PASS(obj) :: MakeZeros => obj_MakeZeros + PROCEDURE, PASS(obj) :: ImportFromToml1 => obj_ImportFromToml1 + PROCEDURE, PASS(obj) :: ImportFromToml2 => obj_ImportFromToml2 + GENERIC, PUBLIC :: ImportFromToml => ImportFromToml1, ImportFromToml2 +END TYPE TDGAlgorithm2_ + +!---------------------------------------------------------------------------- +! IsInitiated@ConstructorMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2025-11-27 +! summary: Check if the object is initiated + +INTERFACE + MODULE FUNCTION obj_IsInitiated(obj) RESULT(ans) + CLASS(TDGAlgorithm2_), INTENT(IN) :: obj + LOGICAL(LGT) :: ans + END FUNCTION obj_IsInitiated +END INTERFACE + +!---------------------------------------------------------------------------- +! Set@ConstructorMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2025-11-06 +! summary: Initiate Newmark-Beta method + +INTERFACE + MODULE SUBROUTINE obj_Initiate(obj, elemsd, facetElemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd, facetElemsd + END SUBROUTINE obj_Initiate +END INTERFACE + +!---------------------------------------------------------------------------- +! MakeZeros@ConstructorMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2025-11-06 +! summary: Reset the TDGAlgorithm2_ object to zero values + +INTERFACE + MODULE SUBROUTINE obj_MakeZeros(obj) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + ! internal varibales + REAL(DFP), PARAMETER :: myzero = 0.0_DFP + END SUBROUTINE obj_MakeZeros +END INTERFACE + +!---------------------------------------------------------------------------- +! Display@ConstructorMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2025-11-06 +! summary: Display the content + +INTERFACE + MODULE SUBROUTINE obj_Display(obj, msg, unitno) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: msg + INTEGER(I4B), OPTIONAL, INTENT(IN) :: unitno + END SUBROUTINE obj_Display +END INTERFACE + +!---------------------------------------------------------------------------- +! Deallocate@ConstructorMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2025-11-06 +! summary: Deallocate the object + +INTERFACE + MODULE SUBROUTINE obj_Deallocate(obj) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + END SUBROUTINE obj_Deallocate +END INTERFACE + +!---------------------------------------------------------------------------- +! ImportFromToml@TomlMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2025-11-06 +! summary: Import data from toml table + +INTERFACE + MODULE SUBROUTINE obj_ImportFromToml1(obj, table) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(toml_table), INTENT(INOUT) :: table + END SUBROUTINE obj_ImportFromToml1 +END INTERFACE + +!---------------------------------------------------------------------------- +! ImportFromToml@TomlMethods +!---------------------------------------------------------------------------- + +!> author: Vikas Sharma, Ph. D. +! date: 2023-11-08 +! summary: Initiate kernel from the toml file + +INTERFACE + MODULE SUBROUTINE obj_ImportFromToml2( & + obj, tomlName, afile, filename, printToml) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + CHARACTER(*), INTENT(IN) :: tomlName + TYPE(TxtFile_), OPTIONAL, INTENT(INOUT) :: afile + CHARACTER(*), OPTIONAL, INTENT(IN) :: filename + LOGICAL(LGT), OPTIONAL, INTENT(IN) :: printToml + ! internal variables + CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml2()" + TYPE(toml_table), ALLOCATABLE :: table + TYPE(toml_table), POINTER :: node + INTEGER(I4B) :: origin, stat + END SUBROUTINE obj_ImportFromToml2 +END INTERFACE + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +END MODULE TDGAlgorithm2_Class diff --git a/src/modules/VectorField/src/VectorField_Class.F90 b/src/modules/VectorField/src/VectorField_Class.F90 index 02228e59..d6adc352 100644 --- a/src/modules/VectorField/src/VectorField_Class.F90 +++ b/src/modules/VectorField/src/VectorField_Class.F90 @@ -164,7 +164,7 @@ MODULE VectorField_Class ! SET: ! @PointNBCMethods - PROCEDURE, NON_OVERRIDABLE, PASS(obj) :: ApplyPointNeumannBC => & + PROCEDURE, PUBLIC, NON_OVERRIDABLE, PASS(obj) :: ApplyPointNeumannBC => & obj_ApplyPointNeumannBC !! Apply point neumann boundary condition diff --git a/src/submodules/STVectorField/CMakeLists.txt b/src/submodules/STVectorField/CMakeLists.txt index 9bdd111f..2606409f 100644 --- a/src/submodules/STVectorField/CMakeLists.txt +++ b/src/submodules/STVectorField/CMakeLists.txt @@ -17,10 +17,15 @@ set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/") target_sources( - ${PROJECT_NAME} - PRIVATE ${src_path}/STVectorField_Class@ConstructorMethods.F90 - ${src_path}/STVectorField_Class@DBCMethods.F90 - ${src_path}/STVectorField_Class@GetMethods.F90 - ${src_path}/STVectorField_Class@IOMethods.F90 - ${src_path}/STVectorField_Class@HDFMethods.F90 - ${src_path}/STVectorField_Class@SetMethods.F90) + ${PROJECT_NAME} + PRIVATE + ${src_path}/STVectorField_Class@BodySourceMethods.F90 + ${src_path}/STVectorField_Class@ConstructorMethods.F90 + ${src_path}/STVectorField_Class@DBCMethods.F90 + ${src_path}/STVectorField_Class@GetMethods.F90 + ${src_path}/STVectorField_Class@HDFMethods.F90 + ${src_path}/STVectorField_Class@IOMethods.F90 + ${src_path}/STVectorField_Class@PointNBCMethods.F90 + ${src_path}/STVectorField_Class@SetMethods.F90 + ${src_path}/STVectorField_Class@SurfaceNBCMethods.F90 +) diff --git a/src/submodules/STVectorField/src/STVectorField_Class@BodySourceMethods.F90 b/src/submodules/STVectorField/src/STVectorField_Class@BodySourceMethods.F90 new file mode 100644 index 00000000..970bd16b --- /dev/null +++ b/src/submodules/STVectorField/src/STVectorField_Class@BodySourceMethods.F90 @@ -0,0 +1,379 @@ +! 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 +! + +SUBMODULE(STVectorField_Class) BodySourceMethods +USE ReallocateUtility, ONLY: Reallocate +USE FEVariable_Method, ONLY: NodalVariable +USE FEVariable_Method, ONLY: QuadratureVariable +USE FEVariable_Method, ONLY: FEVariable_Set => Set +USE FEVariable_Method, ONLY: FEVariable_Deallocate => DEALLOCATE +USE AbstractFE_Class, ONLY: AbstractFE_ +USE AbstractOneDimFE_Class, ONLY: AbstractOneDimFE_ +USE AbstractMesh_Class, ONLY: AbstractMesh_ +USE ForceVector_Method, ONLY: ForceVector_ +USE STForceVector_Method, ONLY: STForceVector_ +USE BaseType, ONLY: QuadraturePoint_ +USE BaseType, ONLY: ElemshapeData_ +USE BaseType, ONLY: FEVariable_ +USE BaseType, ONLY: TypeFEVariableScalar, TypeFEVariableVector +USE BaseType, ONLY: TypeFEVariableSpace +USE BaseType, ONLY: TypeFEVariableSpaceTime +USE BaseType, ONLY: math => TypeMathOpt +USE FieldOpt_Class, ONLY: TypeFieldOpt +USE QuadraturePoint_Method, ONLY: QuadraturePoint_Deallocate => DEALLOCATE +USE ElemshapeData_Method, ONLY: ElemshapeData_Deallocate => DEALLOCATE +USE SwapUtility, ONLY: Swap_ + +#ifdef DEBUG_VER +USE Display_Method, ONLY: Display +USE QuadraturePoint_Method, ONLY: QuadraturePoint_Display => Display +USE ElemshapeData_Method, ONLY: ElemshapeData_Display => Display +USE FEVariable_Method, ONLY: Display +#endif + +IMPLICIT NONE +CONTAINS + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplyBodySource1 +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ApplyBodySource1()" +#endif + +INTEGER(I4B) :: iel, tElements, maxNNS, maxNNT, maxNNSGeo, tcellCon, xij_i, & + xij_j, maxNips, ips, maxNipt, ipt, & + forceVec_i, forceVec_j, forceVec_k, & + spaceCompo(1), itcompo +INTEGER(I4B), ALLOCATABLE :: cellcon(:) +REAL(DFP) :: args(4) +REAL(DFP), ALLOCATABLE :: xij(:, :), forceVec(:, :, :), forceVecQuad(:, :, :), & + forceVec_swap(:, :, :) +TYPE(QuadraturePoint_) :: quad, timeQuad +TYPE(ElemshapeData_) :: elemsd, geoelemsd, timeelemsd, timegeoelemsd +TYPE(FEVariable_) :: forceVar +CLASS(AbstractFE_), POINTER :: feptr, geofeptr +CLASS(AbstractOneDimFE_), POINTER :: timefeptr +CLASS(AbstractMesh_), POINTER :: mesh + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +mesh => obj%fedof%GetMeshPointer() +tElements = mesh%GetTotalElements() + +CALL obj%timefedof%SetFE() +timefeptr => obj%timefedof%GetFEPointer() + +CALL timefeptr%GetGlobalTimeElemShapeData( & + elemsd=timeelemsd, geoelemsd=timegeoelemsd, quad=timeQuad, times=times) + +maxNNSGeo = obj%geofedof%GetMaxTotalConnectivity() +maxNNS = obj%fedof%GetMaxTotalConnectivity() +maxNNT = obj%timefedof%GetMaxTotalConnectivity() +maxNips = obj%fedof%GetMaxTotalQuadraturePoints() +maxNipt = obj%timefedof%GetMaxTotalQuadraturePoints() +spaceCompo = obj%GetspaceCompo(1) + +CALL Reallocate(xij, 3, maxNNSGeo) +CALL Reallocate(cellcon, maxNNS) +CALL Reallocate(forceVec, spaceCompo(1), maxNNS, maxNNT) +CALL Reallocate(forceVec_swap, spaceCompo(1), maxNNT, maxNNS) +CALL Reallocate(forceVecQuad, spaceCompo(1), maxNips, maxNipt) + +forceVar = QuadratureVariable( & + dim1=spaceCompo(1), dim2=maxNips, dim3=maxNipt, & + rank=TypeFEVariableVector, & + vartype=TypeFEVariableSpaceTime) + +args = 0.0_DFP + +DO iel = 1, tElements + + CALL obj%fedof%SetFE(globalElement=iel, islocal=math%yes) + CALL obj%geofedof%SetFE(globalElement=iel, islocal=math%yes) + + feptr => obj%fedof%GetFEPointer(globalElement=iel, islocal=math%yes) + geofeptr => obj%geofedof%GetFEPointer(globalElement=iel, islocal=math%yes) + + CALL obj%fedof%GetConnectivity_( & + globalElement=iel, islocal=math%yes, ans=cellcon, tsize=tcellCon, & + opt="A") + + CALL mesh%GetNodeCoord( & + nodeCoord=xij, nrow=xij_i, ncol=xij_j, islocal=math%yes, & + globalElement=iel) + + CALL feptr%GetGlobalElemShapeData2( & + geofeptr=geofeptr, elemsd=elemsd, geoelemsd=geoelemsd, xij=xij, & + quad=quad) + + DO ipt = 1, timeelemsd%nips + args(4) = timeelemsd%coord(1, ipt) + DO ips = 1, elemsd%nips + args(1:elemsd%nsd) = elemsd%coord(1:elemsd%nsd, ips) + CALL bodySource%GetVectorValue( & + val=forceVecQuad(1:elemsd%nsd, ips, ipt), & + args=args, n=elemsd%nsd) + END DO + END DO + + CALL FEVariable_Set( & + obj=forceVar, & + val=forceVecQuad(1:elemsd%nsd, 1:elemsd%nips, 1:timeelemsd%nips), & + rank=TypeFEVariableVector, varType=TypeFEVariableSpaceTime, & + scale=1.0_DFP, addContribution=math%no) + + CALL STForceVector_( & + testSpace=elemsd, testTime=timeelemsd, c=forceVar, & + crank=TypeFEVariableVector, ans=forceVec, & + dim1=forceVec_i, dim2=forceVec_j, dim3=forceVec_k) + + CALL Swap_(a=forceVec_swap, b=forceVec, i1=1, i2=3, i3=2) + + CALL obj%Set(globalNode=cellcon(1:tcellCon), islocal=math%yes, & + scale=scale, addContribution=math%yes, & + VALUE=forceVec_swap(1:forceVec_i, 1:forceVec_k, 1:forceVec_j)) + +END DO + +DEALLOCATE (cellcon, xij, forceVec, forceVecQuad) +NULLIFY (feptr, geofeptr, timefeptr, mesh) +CALL QuadraturePoint_Deallocate(quad) +CALL QuadraturePoint_Deallocate(timeQuad) +CALL ElemshapeData_Deallocate(elemsd) +CALL ElemshapeData_Deallocate(geoelemsd) +CALL ElemshapeData_Deallocate(timeelemsd) +CALL ElemshapeData_Deallocate(timegeoelemsd) +CALL FEVariable_Deallocate(forceVar) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_ApplyBodySource1 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplyBodySource2 +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ApplyBodySource2()" +#endif + +CALL e%RaiseError(modName//'::'//myName//' - '// & + '[WIP ERROR] :: This routine is under development') + +! LOGICAL(LGT), PARAMETER :: yes = .TRUE., no = .FALSE. +! TYPE(QuadraturePoint_) :: quad +! TYPE(ElemshapeData_) :: elemsd, geoelemsd +! TYPE(FEVariable_) :: forceVar +! INTEGER(I4B) :: iel, tElements, maxNNE, maxNNEGeo, & +! tcellCon, tforceVec, xij_i, xij_j, & +! tfevec +! INTEGER(I4B), ALLOCATABLE :: cellcon(:) +! REAL(DFP), ALLOCATABLE :: xij(:, :), fevec(:), forceVec(:) +! CLASS(AbstractFE_), POINTER :: feptr, geofeptr +! CLASS(AbstractMesh_), POINTER :: mesh +! +! #ifdef DEBUG_VER +! CALL e%RaiseInformation(modName//'::'//myName//' - '// & +! '[START] ') +! #endif +! +! mesh => obj%fedof%GetMeshPointer() +! tElements = mesh%GetTotalElements() +! maxNNEGeo = obj%geofedof%GetMaxTotalConnectivity() +! maxNNE = obj%fedof%GetMaxTotalConnectivity() +! +! CALL Reallocate(xij, 3, maxNNEGeo) +! CALL Reallocate(cellcon, maxNNE) +! CALL Reallocate(fevec, maxNNE) +! CALL Reallocate(forceVec, maxNNE) +! +! forceVar = NodalVariable(tsize=maxNNE, rank=TypeFEVariableScalar, & +! vartype=TypeFEVariableSpace) +! ! forceVar is nodal value at element level obtained from the bodySource +! ! field +! ! TODO: +! ! We should get the quadrature values of bodySource in forceVar. +! ! This means we should add a method to scalarfield for getting the +! ! quadrature values of a field. This means, the method will first create +! ! local element shape data, and then interpolate the data using nodal +! ! values and the shape data. This interpolated value should be returned. +! ! In this way, we can use any order of interpolation for the bodySource. +! +! DO iel = 1, tElements +! +! CALL obj%fedof%SetFE(globalElement=iel, islocal=yes) +! CALL obj%geofedof%SetFE(globalElement=iel, islocal=yes) +! +! feptr => obj%fedof%GetFEPointer(globalElement=iel, islocal=yes) +! geofeptr => obj%geofedof%GetFEPointer(globalElement=iel, islocal=yes) +! +! CALL obj%fedof%GetConnectivity_(globalElement=iel, islocal=yes, & +! ans=cellcon, tsize=tcellCon, opt="A") +! +! CALL mesh%GetNodeCoord(nodeCoord=xij, nrow=xij_i, ncol=xij_j, & +! islocal=yes, globalElement=iel) +! +! CALL feptr%GetGlobalElemShapeData2( & +! geofeptr=geofeptr, elemsd=elemsd, geoelemsd=geoelemsd, xij=xij, & +! quad=quad) +! +! ! Read the above TODO comment +! CALL bodySource%Get(VALUE=forceVec, globalNode=cellcon(1:tcellCon), & +! tsize=tforceVec, islocal=yes) +! +! CALL FEVariable_Set( & +! obj=forceVar, val=forceVec(1:tforceVec), rank=TypeFEVariableScalar, & +! vartype=TypeFEVariableSpace, scale=1.0_DFP, addContribution=no) +! +! CALL ForceVector_( & +! test=elemsd, c=forceVar, crank=TypeFEVariableScalar, ans=fevec, & +! tsize=tfevec) +! +! CALL obj%Set(globalNode=cellcon(1:tcellCon), islocal=yes, & +! scale=scale, addContribution=yes, & +! VALUE=fevec(1:tfevec)) +! +! END DO +! +! IF (ALLOCATED(cellcon)) DEALLOCATE (cellcon) +! IF (ALLOCATED(xij)) DEALLOCATE (xij) +! IF (ALLOCATED(forceVec)) DEALLOCATE (forceVec) +! IF (ALLOCATED(fevec)) DEALLOCATE (fevec) +! NULLIFY (feptr, geofeptr, mesh) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_ApplyBodySource2 + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplyBodySource3 +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ApplyBodySource3()" +#endif + +LOGICAL(LGT), PARAMETER :: yes = .TRUE., no = .FALSE. + +INTEGER(I4B) :: iel, tElements, maxNNE, maxNNEGeo, & + tcellCon, tforceVec, xij_i, xij_j, maxNips, ips, & + spaceCompo(1), force_i, force_j +REAL(DFP) :: args(4) + +TYPE(QuadraturePoint_) :: quad +TYPE(ElemshapeData_) :: elemsd, geoelemsd +TYPE(FEVariable_) :: forceVar +CLASS(AbstractFE_), POINTER :: feptr, geofeptr +CLASS(AbstractMesh_), POINTER :: mesh + +REAL(DFP), ALLOCATABLE :: xij(:, :), forceVec(:, :), forceVecQuad(:, :) +INTEGER(I4B), ALLOCATABLE :: cellcon(:) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +mesh => obj%fedof%GetMeshPointer() + +tElements = mesh%GetTotalElements() +maxNNEGeo = obj%geofedof%GetMaxTotalConnectivity() +maxNNE = obj%fedof%GetMaxTotalConnectivity() +maxNips = obj%fedof%GetMaxTotalQuadraturePoints() +spaceCompo = obj%GetspaceCompo(1) + +CALL Reallocate(xij, 3, maxNNEGeo) +CALL Reallocate(cellcon, maxNNE) +CALL Reallocate(forceVec, spaceCompo(1), maxNNE) +CALL Reallocate(forceVecQuad, spaceCompo(1), maxNips) + +forceVar = QuadratureVariable(nrow=spaceCompo(1), & + ncol=maxNips, rank=TypeFEVariableVector, & + vartype=TypeFEVariableSpace) + +args = 0.0_DFP +args(4) = times + +DO iel = 1, tElements + + CALL obj%fedof%SetFE(globalElement=iel, islocal=yes) + CALL obj%geofedof%SetFE(globalElement=iel, islocal=yes) + + feptr => obj%fedof%GetFEPointer(globalElement=iel, islocal=yes) + geofeptr => obj%geofedof%GetFEPointer(globalElement=iel, islocal=yes) + + CALL obj%fedof%GetConnectivity_(globalElement=iel, islocal=yes, & + ans=cellcon, tsize=tcellCon, opt="A") + + CALL mesh%GetNodeCoord( & + nodeCoord=xij, nrow=xij_i, ncol=xij_j, islocal=yes, globalElement=iel) + + CALL feptr%GetGlobalElemShapeData2( & + geofeptr=geofeptr, elemsd=elemsd, geoelemsd=geoelemsd, xij=xij, & + quad=quad) + + DO ips = 1, elemsd%nips + args(1:elemsd%nsd) = elemsd%coord(1:elemsd%nsd, ips) + CALL bodySource%GetVectorValue(val=forceVecQuad(1:elemsd%nsd, ips), & + args=args, n=elemsd%nsd) + END DO + + CALL FEVariable_Set( & + obj=forceVar, val=forceVecQuad(1:elemsd%nsd, 1:elemsd%nips), & + rank=TypeFEVariableVector, varType=TypeFEVariableSpace, & + scale=1.0_DFP, addContribution=no) + + CALL ForceVector_(test=elemsd, c=forceVar, crank=TypeFEVariableVector, & + ans=forceVec, nrow=force_i, ncol=force_j) + + CALL obj%Set( & + globalNode=cellcon(1:tcellCon), islocal=yes, scale=scale, & + addContribution=yes, VALUE=forceVec(1:force_i, 1:force_j)) + +END DO + +IF (ALLOCATED(xij)) DEALLOCATE (xij) +IF (ALLOCATED(cellcon)) DEALLOCATE (cellcon) +IF (ALLOCATED(forceVec)) DEALLOCATE (forceVec) +IF (ALLOCATED(forceVecQuad)) DEALLOCATE (forceVecQuad) +NULLIFY (feptr, geofeptr, mesh) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_ApplyBodySource3 + +!---------------------------------------------------------------------------- +! Include Errors +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE BodySourceMethods diff --git a/src/submodules/STVectorField/src/STVectorField_Class@ConstructorMethods.F90 b/src/submodules/STVectorField/src/STVectorField_Class@ConstructorMethods.F90 index 9e11d20f..2a772ece 100644 --- a/src/submodules/STVectorField/src/STVectorField_Class@ConstructorMethods.F90 +++ b/src/submodules/STVectorField/src/STVectorField_Class@ConstructorMethods.F90 @@ -65,6 +65,82 @@ END SELECT END PROCEDURE obj_Initiate2 +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Initiate4 +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Initiate4()" +#endif + +CHARACTER(1) :: dof_names(1) +INTEGER(I4B) :: dof_tNodes(1), dof_tsize, dof_spaceCompo(1), & + dof_timeCompo(1), dof_tPhysicalVarNames +LOGICAL(LGT) :: istimefedof, isok, timeCompoMade + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +CALL obj%DEALLOCATE() + +timeCompoMade = .FALSE. +istimefedof = PRESENT(timefedof) +IF (istimefedof) THEN + timeCompoMade = timefedof%IsInitiated() + IF (timeCompoMade) dof_timeCompo(1) = timefedof%GetTotalDOF() +END IF + +IF (.NOT. timeCompoMade) THEN +#ifdef DEBUG_VER + isok = PRESENT(timeCompo) + CALL AssertError1(isok, myName, "timeCompo is not present") +#endif + dof_timeCompo(1) = timeCompo(1) + timeCompoMade = .TRUE. +END IF + +isok = PRESENT(spaceCompo) +CALL AssertError1(isok, myName, "spaceCompo is not present") + +dof_spaceCompo(1) = spaceCompo(1) + +obj%timeCompo = dof_timeCompo(1) +obj%spaceCompo = dof_spaceCompo(1) + +CALL Reallocate(obj%space_idofs, obj%spaceCompo) +obj%space_idofs = Arange(1_I4B, obj%spaceCompo) + +CALL Reallocate(obj%time_idofs, obj%timeCompo) +obj%time_idofs = Arange(1_I4B, obj%timeCompo) + +CALL Reallocate(obj%idofs, obj%timeCompo * obj%spaceCompo) +obj%idofs = Arange(1_I4B, obj%timeCompo * obj%spaceCompo) + +dof_names(1) = name(1:1) +dof_tNodes(1) = fedof%GetTotalDOF() +dof_tsize = dof_tNodes(1) * dof_spaceCompo(1) * dof_timeCompo(1) +dof_tPhysicalVarNames = 1_I4B + +CALL AbstractNodeFieldInitiate( & + obj=obj, name=name, engine=engine, fieldType=fieldType, comm=comm, & + local_n=local_n, global_n=global_n, fedof=fedof, timefedof=timefedof, & + storageFMT=MYSTORAGEFORMAT, spaceCompo=dof_spaceCompo, & + isSpaceCompo=.TRUE., isSpaceCompoScalar=.TRUE., timeCompo=dof_timeCompo, & + isTimeCompo=.TRUE., isTimeCompoScalar=.TRUE., & + tPhysicalVarNames=dof_tPhysicalVarNames, physicalVarNames=dof_names, & + isPhysicalVarNames=.TRUE., isPhysicalVarNamesScalar=.TRUE., & + tSize=dof_tsize, tNodes=dof_tNodes, isTNodes=.TRUE., & + isTNodesScalar=.TRUE., geofedof=geofedof) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_Initiate4 + !---------------------------------------------------------------------------- ! Deallocate !---------------------------------------------------------------------------- diff --git a/src/submodules/STVectorField/src/STVectorField_Class@PointNBCMethods.F90 b/src/submodules/STVectorField/src/STVectorField_Class@PointNBCMethods.F90 new file mode 100644 index 00000000..fd623f8a --- /dev/null +++ b/src/submodules/STVectorField/src/STVectorField_Class@PointNBCMethods.F90 @@ -0,0 +1,102 @@ +! 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 +! + +SUBMODULE(STVectorField_Class) PointNBCMethods +USE Display_Method, ONLY: ToString +USE ReallocateUtility, ONLY: Reallocate +USE NeumannBC_Class, ONLY: NeumannBC_ +USE BaseType, ONLY: math => TypeMathOpt + +#ifdef DEBUG_VER +USE Display_Method, ONLY: Display +#endif + +IMPLICIT NONE +CONTAINS + +!---------------------------------------------------------------------------- +! ApplyPointNeumannBC +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplyPointNeumannBC +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ApplyPointNeumannBC1()" +#endif + +INTEGER(I4B), PARAMETER :: expandFactor = 2 +INTEGER(I4B) :: idof, nrow, ncol, tbc, ibc, spaceCompo +CLASS(NeumannBC_), POINTER :: nbcptr +LOGICAL(LGT) :: isok + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +ncol = obj%timeCompo +nrow = obj%GetMaxTotalNodeNumForBC() + +CALL Reallocate(obj%nodalvalue, nrow, ncol, isExpand=math%yes, & + expandFactor=expandFactor) +CALL Reallocate(obj%nodenum, nrow, isExpand=math%yes, & + expandFactor=expandFactor) + +tbc = SIZE(obj%nbc_point) +DO ibc = 1, tbc + nbcptr => obj%nbc_point(ibc)%ptr + isok = ASSOCIATED(nbcptr) + IF (.NOT. isok) CYCLE + + spaceCompo = nbcptr%GetDOFNo() + +#ifdef DEBUG_VER + isok = spaceCompo .LE. obj%spaceCompo + CALL AssertError1( & + isok, myName, & + "DOFNo obtained from nbc_point("//ToString(ibc)// & + ") is "//ToString(spaceCompo)//" is greater than total space & + &components of the field which is "//ToString(obj%spaceCompo)) +#endif + + CALL nbcptr%Get( & + nodalvalue=obj%nodalvalue, nodenum=obj%nodenum, times=times, & + nrow=nrow, ncol=ncol, fedof=obj%fedof, geofedof=obj%geofedof, & + timefedof=obj%timefedof) + + DO idof = 1, ncol + CALL obj%Set( & + globalNode=obj%nodenum(1:nrow), VALUE=obj%nodalvalue(1:nrow, idof), & + scale=scale, addContribution=math%yes, islocal=math%yes, & + spaceCompo=spaceCompo, timeCompo=idof) + END DO +END DO + +nbcptr => NULL() + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_ApplyPointNeumannBC + +!---------------------------------------------------------------------------- +! Include Errors +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE PointNBCMethods diff --git a/src/submodules/STVectorField/src/STVectorField_Class@SetMethods.F90 b/src/submodules/STVectorField/src/STVectorField_Class@SetMethods.F90 index 33671b38..2a68c4b0 100644 --- a/src/submodules/STVectorField/src/STVectorField_Class@SetMethods.F90 +++ b/src/submodules/STVectorField/src/STVectorField_Class@SetMethods.F90 @@ -292,13 +292,11 @@ #endif -!$OMP PARALLEL DO PRIVATE(ii) DO ii = 1, SIZE(globalNode) CALL obj%Set(VALUE=VALUE(:, :, ii), globalNode=globalNode(ii), & scale=scale, addContribution=addContribution, & islocal=islocal) END DO -!$OMP END PARALLEL DO END PROCEDURE obj_Set7 !---------------------------------------------------------------------------- @@ -484,7 +482,6 @@ indx = obj%dof.tNodes.1 -!$OMP PARALLEL DO PRIVATE(ii) DO ii = 1, indx CALL GetNodeLoc_(obj=obj%dof, nodenum=ii, ans=TEMP_INTVEC, tsize=tsize, & @@ -494,7 +491,6 @@ scale=scale, addContribution=addContribution) END DO -!$OMP END PARALLEL DO END PROCEDURE obj_Set12 @@ -755,16 +751,92 @@ !---------------------------------------------------------------------------- MODULE PROCEDURE obj_SetFromVectorField -INTEGER(I4B) :: ii, jj +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_SetFromVectorField()" +LOGICAL(LGT) :: isok +#endif + +INTEGER(I4B) :: idof, icompo + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif -DO ii = 1, obj%spaceCompo - jj = GetIDOF(spaceCompo=ii, timeCompo=1, tspaceCompo=obj%spaceCompo) - CALL obj%Set(idof=jj, ivar=1, VALUE=VALUE, idof_value=ii, & +#ifdef DEBUG_VER +isok = obj%IsInitiated() +CALL AssertError1(isok, myName, & + 'STVectorField_::obj is not initiated') +#endif + +#ifdef DEBUG_VER +isok = VALUE%IsInitiated() +CALL AssertError1(isok, myName, & + 'VectorField_::value is not initiated') +#endif + +DO icompo = 1, obj%spaceCompo + idof = GetIDOF(spaceCompo=icompo, timeCompo=timeCompo, tspaceCompo=obj%spaceCompo) + CALL obj%Set(idof=idof, ivar=1, VALUE=VALUE, idof_value=icompo, & ivar_value=1, scale=scale, addContribution=addContribution) END DO END PROCEDURE obj_SetFromVectorField +!---------------------------------------------------------------------------- +! SetToVectorField +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_SetToVectorField +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_SetToVectorField()" +LOGICAL(LGT) :: isok +#endif + +INTEGER(I4B) :: s(3), p(3), icompo, idof +REAL(DFP), POINTER :: realvec(:) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +#ifdef DEBUG_VER +isok = obj%IsInitiated() +CALL AssertError1(isok, myName, & + 'STVectorField_::obj is not initiated') +#endif + +#ifdef DEBUG_VER +isok = VALUE%IsInitiated() +CALL AssertError1(isok, myName, & + 'VectorField_::value is not initiated') +#endif + +realvec => obj%GetPointer() + +#ifdef DEBUG_VER +isok = ASSOCIATED(realvec) +CALL AssertError1(isok, myName, & + 'realvec obtained To value is not ASSOCIATED') +#endif + +DO icompo = 1, obj%spaceCompo + s = GetNodeLoc(obj=VALUE%dof, idof=icompo) + + idof = GetIDOF(spaceCompo=icompo, timeCompo=timeCompo, tspaceCompo=obj%spaceCompo) + p = GetNodeLoc(obj=obj%dof, idof=idof) + + CALL VALUE%SetMultiple( & + VALUE=realvec, scale=scale, addContribution=addContribution, & + istart_value=p(1), iend_value=p(2), stride_value=p(3), & + istart=s(1), iend=s(2), stride=s(3)) +END DO + +realvec => NULL() + +END PROCEDURE obj_SetToVectorField + !---------------------------------------------------------------------------- ! !---------------------------------------------------------------------------- diff --git a/src/submodules/STVectorField/src/STVectorField_Class@SurfaceNBCMethods.F90 b/src/submodules/STVectorField/src/STVectorField_Class@SurfaceNBCMethods.F90 new file mode 100644 index 00000000..eea7ca43 --- /dev/null +++ b/src/submodules/STVectorField/src/STVectorField_Class@SurfaceNBCMethods.F90 @@ -0,0 +1,292 @@ +! 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 +! + +SUBMODULE(STVectorField_Class) SurfaceNBCMethods +USE BaseType, ONLY: QuadraturePoint_ +USE BaseType, ONLY: ElemShapeData_ +USE BaseType, ONLY: TypeFEVariableVector +USE BaseType, ONLY: TypeFEVariableSpaceTime +USE BaseType, ONLY: math => TypeMathOpt +USE ReallocateUtility, ONLY: Reallocate +USE FEVariable_Method, ONLY: NodalVariable +USE FEVariable_Method, ONLY: FEVariable_Set => Set +USE FEVariable_Method, ONLY: FEVariable_Deallocate => DEALLOCATE +USE QuadraturePoint_Method, ONLY: QuadraturePoint_Deallocate => DEALLOCATE +USE ElemshapeData_Method, ONLY: ElemShapeData_Deallocate => DEALLOCATE +USE AbstractFE_Class, ONLY: AbstractFE_ +USE AbstractOneDimFE_Class, ONLY: AbstractOneDimFE_ +USE STForceVector_Method, ONLY: STForceVector_ +USE NeumannBC_Class, ONLY: NeumannBC_ +USE AbstractMesh_Class, ONLY: AbstractMesh_ +USE FieldOpt_Class, ONLY: TypeFieldOpt +USE SwapUtility, ONLY: Swap_ + +#ifdef DEBUG_VER +USE Display_Method, ONLY: Display +USE QuadraturePoint_Method, ONLY: QuadraturePoint_Display => Display +USE ElemshapeData_Method, ONLY: ElemshapeData_Display => Display +#endif + +IMPLICIT NONE +CONTAINS + +!---------------------------------------------------------------------------- +! ApplySurfaceNeumannBC +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ApplySurfaceNeumannBC +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ApplySurfaceNeumannBC()" +#endif + +LOGICAL(LGT) :: isok + +CLASS(NeumannBC_), POINTER :: nbc +CLASS(AbstractMesh_), POINTER :: mesh +CLASS(AbstractOneDimFE_), POINTER :: timefeptr + +INTEGER(I4B) :: tbc, ibc, maxNNSGeo, maxNNS, maxNNT, spaceCompo(1) +INTEGER(I4B), ALLOCATABLE :: facetCon(:) +REAL(DFP), ALLOCATABLE :: xij(:, :), nbcValue(:, :, :), forceVec(:, :, :) + +TYPE(FEVariable_) :: forceVar +TYPE(QuadraturePoint_) :: quad, facetQuad, timeQuad +TYPE(ElemShapeData_) :: elemsd, facetElemsd, geoElemsd, geoFacetElemsd, & + timeelemsd, timegeoelemsd + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +tbc = obj%GetTotalNBC() +isok = tbc .GT. 0 + +IF (.NOT. isok) THEN +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif + RETURN +END IF + +CALL nbcField%SetAll(VALUE=math%zero) + +DO ibc = 1, tbc + nbc => obj%GetNBCPointer(ibc) + isok = ASSOCIATED(nbc) + IF (.NOT. isok) CYCLE + CALL nbcField%ApplyDirichletBC(dbc=nbc, times=times) +END DO +nbc => NULL() + +mesh => obj%fedof%GetMeshPointer() +CALL obj%timefedof%SetFE() +timefeptr => obj%timefedof%GetFEPointer() +CALL timefeptr%GetGlobalTimeElemShapeData( & + elemsd=timeelemsd, geoelemsd=timegeoelemsd, quad=timeQuad, times=times) + +maxNNSGeo = obj%geofedof%GetMaxTotalConnectivity() +maxNNS = obj%fedof%GetMaxTotalConnectivity() +maxNNT = obj%timefedof%GetMaxTotalConnectivity() +spaceCompo = obj%GetSpaceCompo(1) + +CALL Reallocate(xij, 3, maxNNSGeo) +CALL Reallocate(facetCon, maxNNS) +CALL Reallocate(forceVec, spaceCompo(1), maxNNS, maxNNT) +CALL Reallocate(nbcValue, spaceCompo(1), maxNNT, maxNNS) ! nodes FMT + +forceVar = NodalVariable( & + dim1=spaceCompo(1), dim2=maxNNT, dim3=maxNNS, & + rank=TypeFEVariableVector, & + vartype=TypeFEVariableSpaceTime) + +DO ibc = 1, tbc + nbc => obj%GetNBCPointer(ibc) + + isok = ASSOCIATED(nbc) + IF (.NOT. isok) CYCLE + + CALL STVectorFieldAssembleSurfaceSource( & + obj=obj, nbc=nbc, fedof=obj%fedof, geofedof=obj%geofedof, & + mesh=mesh, nbcField=nbcField, scale=scale, xij=xij, & + forceVec=forceVec, nbcValue=nbcValue, forceVar=forceVar, quad=quad, & + facetQuad=facetQuad, elemsd=elemsd, facetElemsd=facetElemsd, & + geoElemsd=geoElemsd, geoFacetElemsd=geoFacetElemsd, & + facetCon=facetCon, timeelemsd=timeelemsd) +END DO + +DEALLOCATE (xij, nbcValue, forceVec, facetCon) +NULLIFY (timefeptr, mesh, nbc) +CALL FEVariable_Deallocate(forceVar) +CALL QuadraturePoint_Deallocate(quad) +CALL QuadraturePoint_Deallocate(facetQuad) +CALL QuadraturePoint_Deallocate(timeQuad) +CALL ElemshapeData_Deallocate(elemsd) +CALL ElemshapeData_Deallocate(facetElemsd) +CALL ElemshapeData_Deallocate(geoElemsd) +CALL ElemshapeData_Deallocate(geoFacetElemsd) +CALL ElemshapeData_Deallocate(timeelemsd) +CALL ElemshapeData_Deallocate(timegeoelemsd) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_ApplySurfaceNeumannBC + +!---------------------------------------------------------------------------- +! STScalarFieldAssembleSurfaceSource@ScalarFieldMethods +!---------------------------------------------------------------------------- + +!> author: Shion Shimizu +! date: 2025-12-11 +! summary: Assemble surface source +! +!# Introduction +! +! This is an internal routine that assemble surface source (Neumann BC) +! into the right-hand side (RHS) vector of a vector field. + +SUBROUTINE STVectorFieldAssembleSurfaceSource( & + obj, nbc, fedof, geofedof, mesh, nbcField, scale, xij, forceVec, & + nbcValue, forceVar, quad, facetQuad, elemsd, facetElemsd, geoElemsd, & + geoFacetElemsd, facetCon, timeelemsd) + CLASS(STVectorField_), INTENT(INOUT) :: obj + CLASS(NeumannBC_), INTENT(INOUT) :: nbc + CLASS(FEDOF_), INTENT(INOUT) :: fedof + CLASS(FEDOF_), INTENT(INOUT) :: geofedof + CLASS(AbstractMesh_), INTENT(INOUT) :: mesh + CLASS(STVectorField_), INTENT(INOUT) :: nbcField + REAL(DFP), INTENT(IN) :: scale + INTEGER(I4B), INTENT(INOUT) :: facetCon(:) + !! Working arrays for connectivity + !! geoCellCon: connectivity of the geometry + !! geoFacetCon: connectivity of the geometry facet element + !! cellCon: connectivity of the call + REAL(DFP), INTENT(INOUT) :: xij(:, :), forceVec(:, :, :), nbcValue(:, :, :) + !! Working arrays for element vectors and matrix + TYPE(FEVariable_), INTENT(INOUT) :: forceVar + !! Working variables + TYPE(QuadraturePoint_), INTENT(INOUT) :: quad, facetQuad + !! Working variables + TYPE(ElemShapeData_), INTENT(INOUT) :: elemsd, facetElemsd, geoElemsd, & + geoFacetElemsd, timeelemsd + !! Working variables for containing shape function data + +#ifdef DEBUG_VER + CHARACTER(*), PARAMETER :: myName = "STVectorFieldAssembleSurfaceSource()" +#endif + + REAL(DFP) :: forceVec_swap(SIZE(forceVec, 1), SIZE(forceVec, 3), & + SIZE(forceVec, 2)) + LOGICAL(LGT) :: isElemToEdge, isElemToFace, isok + INTEGER(I4B) :: tElemToFace, indx, localCellNumber, localFaceNumber, & + forceVec_i, forceVec_j, forceVec_k, & + xij_i, xij_j, tFacetCon, & + nbcValue_i, nbcValue_j, nbcValue_k + CLASS(AbstractFE_), POINTER :: feptr, geofeptr + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + + isElemToEdge = nbc%IsElemToEdgeInitiated() + isElemToFace = nbc%IsElemToFaceInitiated() + + isok = isElemToEdge .OR. isElemToFace + IF (.NOT. isok) THEN + CALL e%RaiseDebug(modName//'::'//myName//' - '// & + 'isElemToEdge and isElemToFace are both .false. So, nothing to do.') + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif + + RETURN + END IF + + tElemToFace = nbc%GetTotalElemToFace() + + DO indx = 1, tElemToFace + CALL nbc%GetElemToFace(indx=indx, localCellNumber=localCellNumber, & + localFaceNumber=localFaceNumber) + + CALL fedof%SetFE(globalElement=localCellNumber, islocal=math%yes) + feptr => fedof%GetFEPointer(globalElement=localCellNumber, & + islocal=math%yes) + + CALL geofedof%SetFE(globalElement=localCellNumber, islocal=math%yes) + geofeptr => geofedof%GetFEPointer(globalElement=localCellNumber, & + islocal=math%yes) + + CALL fedof%GetFacetConnectivity_( & + globalElement=localCellNumber, islocal=math%yes, ans=facetCon, & + tsize=tFacetCon, localFaceNumber=localFaceNumber) + + CALL mesh%GetNodeCoord( & + nodeCoord=xij, nrow=xij_i, ncol=xij_j, islocal=math%yes, & + globalElement=localCellNumber) + + CALL feptr%GetGlobalFacetElemShapeData2( & + geofeptr=geofeptr, elemsd=elemsd, facetElemsd=facetElemsd, & + geoElemsd=geoElemsd, geoFacetElemsd=geoFacetElemsd, & + localFaceNumber=localFaceNumber, quad=quad, facetQuad=facetQuad, & + xij=xij) + + CALL nbcField%Get( & + VALUE=nbcValue, globalNode=facetCon(1:tFacetCon), & + dim1=nbcValue_i, dim2=nbcValue_j, dim3=nbcValue_k, & + islocal=math%yes, storageFMT=TypeFieldOpt%storageFormatNodes) + + CALL FEVariable_Set( & + obj=forceVar, val=nbcValue(1:nbcValue_i, 1:nbcValue_j, 1:nbcValue_k), & + rank=TypeFEVariableVector, vartype=TypeFEVariableSpaceTime, & + scale=math%one, addContribution=math%no) + + CALL STForceVector_( & + testSpace=facetElemsd, testTime=timeelemsd, c=forceVar, & + crank=TypeFEVariableVector, ans=forceVec, & + dim1=forceVec_i, dim2=forceVec_j, dim3=forceVec_k) + + CALL Swap_(a=forceVec_swap, b=forceVec, i1=1, i2=3, i3=2) + + CALL obj%Set( & + globalNode=facetCon(1:tFacetCon), & + VALUE=forceVec_swap(1:forceVec_i, 1:forceVec_k, 1:forceVec_j), & + scale=scale, addContribution=math%yes, islocal=math%yes) + + END DO + + STOP + NULLIFY (feptr, geofeptr) + +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END SUBROUTINE STVectorFieldAssembleSurfaceSource + +!---------------------------------------------------------------------------- +! Include Errors +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE SurfaceNBCMethods diff --git a/src/submodules/TDGAlgorithms/CMakeLists.txt b/src/submodules/TDGAlgorithms/CMakeLists.txt index 896fe785..0e3a4ec0 100644 --- a/src/submodules/TDGAlgorithms/CMakeLists.txt +++ b/src/submodules/TDGAlgorithms/CMakeLists.txt @@ -17,7 +17,12 @@ set(src_path "${CMAKE_CURRENT_LIST_DIR}/src/") target_sources( - ${PROJECT_NAME} - PRIVATE ${src_path}/TDGAlgorithm1_Class@ConstructorMethods.F90 - ${src_path}/TDGAlgorithm1_Class@IOMethods.F90 - ${src_path}/TDGAlgorithm1_Class@TomlMethods.F90) + ${PROJECT_NAME} + PRIVATE + ${src_path}/TDGAlgorithm1_Class@ConstructorMethods.F90 + ${src_path}/TDGAlgorithm1_Class@IOMethods.F90 + ${src_path}/TDGAlgorithm1_Class@TomlMethods.F90 + ${src_path}/TDGAlgorithm2_Class@ConstructorMethods.F90 + ${src_path}/TDGAlgorithm2_Class@IOMethods.F90 + ${src_path}/TDGAlgorithm2_Class@TomlMethods.F90 +) diff --git a/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@ConstructorMethods.F90 b/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@ConstructorMethods.F90 new file mode 100644 index 00000000..958a3dbe --- /dev/null +++ b/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@ConstructorMethods.F90 @@ -0,0 +1,357 @@ +! 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 +! + +SUBMODULE(TDGAlgorithm2_Class) ConstructorMethods +USE ApproxUtility, ONLY: OPERATOR(.approxeq.) +USE InputUtility, ONLY: Input +USE MassMatrix_Method, ONLY: MassMatrix_ +USE ProductUtility, ONLY: OuterProd_ +USE BaseType, ONLY: math => TypeMathOpt +USE Lapack_Method, ONLY: GetInvMat + +IMPLICIT NONE +CONTAINS + +!---------------------------------------------------------------------------- +! IsInitiated +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_IsInitiated +ans = obj%isInit +END PROCEDURE obj_IsInitiated + +!---------------------------------------------------------------------------- +! Initiate +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Initiate +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Initiate()" +#endif + +INTEGER(I4B) :: i1, i2 + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +CALL obj%DEALLOCATE() + +obj%isInit = .TRUE. + +obj%name(1:4) = "TDG " +obj%nrow = elemsd%nns +obj%ncol = obj%nrow + +CALL GetCt(obj, elemsd) + +CALL GetMt(obj, elemsd, facetElemsd) + +CALL GetWt(obj, elemsd) + +CALL GetAt(obj, elemsd, facetElemsd) + +CALL GetBt(obj, elemsd, facetElemsd) + +CALL GetKt(obj, elemsd) + +obj%dis(1) = obj%at_right + +DO i1 = 1, obj%nrow + obj%dis(i1 + 3) = obj%bt_right(i1) + obj%vel(i1 + 3) = facetElemsd%N(i1, 2) + obj%acc(i1 + 3) = facetElemsd%dNdXt(i1, 1, 2) + obj%rhs_m_v1(i1) = facetElemsd%N(i1, 1) + obj%rhs_k_u1(i1) = -obj%tat(i1) ! minus sign +END DO + +CALL obj%MakeZeros() + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_Initiate + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +! Getting Kt, coefficients for the stiffness matrix +! This should be called after GetBt + +SUBROUTINE GetKt(obj, elemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd + INTEGER(I4B) :: nrow, ncol, ips, ii, jj + REAL(DFP) :: scale + + nrow = obj%nrow + ncol = nrow + obj%kt(1:nrow, 1:ncol) = math%zero + + DO ips = 1, elemsd%nips + + scale = elemsd%ws(ips) * elemsd%thickness(ips) * elemsd%js(ips) + + CALL OuterProd_(a=elemsd%N(1:nrow, ips), & + b=obj%bt(1:ncol, ips), & + ans=obj%kt, nrow=ii, ncol=jj, & + anscoeff=math%one, scale=scale) + + END DO + +END SUBROUTINE GetKt + +!---------------------------------------------------------------------------- +! +!---------------------------------------------------------------------------- + +! Getting Ct, coefficients for the damping matrix + +SUBROUTINE GetCt(obj, elemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd + + CALL MassMatrix_(test=elemsd, trial=elemsd, ans=obj%ct, nrow=obj%nrow, & + ncol=obj%ncol) + +END SUBROUTINE GetCt + +!---------------------------------------------------------------------------- +! GetMt +!---------------------------------------------------------------------------- + +! Getting Mt, coefficients for the mass matrix (not temporal mass matrix) + +SUBROUTINE GetMt(obj, elemsd, facetElemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd, facetElemsd + + INTEGER(I4B) :: i1, i2 + + CALL MassMatrix_(N=elemsd%N, M=elemsd%dNdXt(:, 1, :), & + js=elemsd%js, ws=elemsd%ws, thickness=elemsd%thickness, & + nips=elemsd%nips, nns1=elemsd%nns, nns2=elemsd%nns, & + ans=obj%mt, nrow=obj%nrow, ncol=obj%ncol) + + CALL OuterProd_(a=facetElemsd%N(1:obj%nrow, 1), & + b=facetElemsd%N(1:obj%nrow, 1), & + ans=obj%mt, nrow=i1, ncol=i2, scale=math%one, & + anscoeff=math%one) + +END SUBROUTINE GetMt + +!---------------------------------------------------------------------------- +! GetWt +!---------------------------------------------------------------------------- + +! Getting Wt and Wmt. This should be called after GetMt + +SUBROUTINE GetWt(obj, elemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd + + INTEGER(I4B) :: nrow, ncol + + nrow = obj%nrow + ncol = obj%ncol + + obj%wt(1:nrow, 1:ncol) = obj%mt(1:nrow, 1:ncol) ! temp mass + jump contri + CALL GetInvMat(obj%wt(1:nrow, 1:ncol)) + obj%wmt(1:nrow, 1:ncol) = & + MATMUL(obj%wt(1:nrow, 1:ncol), obj%ct(1:nrow, 1:ncol)) + obj%wmt(1:nrow, 1:ncol) = TRANSPOSE(obj%wmt(1:nrow, 1:ncol)) + +END SUBROUTINE GetWt + +!---------------------------------------------------------------------------- +! GetAt +!---------------------------------------------------------------------------- + +! Getting At +! This should be called after GetWt + +SUBROUTINE GetAt(obj, elemsd, facetElemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd, facetElemsd + REAL(DFP) :: temp(MAX_ORDER_TIME), areal + INTEGER(I4B) :: nrow, ncol, tsize, ii + + nrow = obj%nrow + ncol = obj%ncol + + tsize = elemsd%nips + + temp(1:nrow) = & + MATMUL(obj%wt(1:nrow, 1:ncol), facetElemsd%N(1:ncol, 1)) + + obj%tat(1:nrow) = math%zero + + obj%at_right = DOT_PRODUCT(facetElemsd%N(1:nrow, 2), temp(1:nrow)) + + DO ii = 1, tsize + obj%at(ii) = DOT_PRODUCT(elemsd%N(1:nrow, ii), temp(1:nrow)) + areal = obj%at(ii) * elemsd%ws(ii) * elemsd%thickness(ii) * elemsd%js(ii) + obj%tat(1:nrow) = obj%tat(1:nrow) + elemsd%N(1:nrow, ii) * areal + END DO + +END SUBROUTINE GetAt + +!---------------------------------------------------------------------------- +! GetBt +!---------------------------------------------------------------------------- + +! Getting Bt +! This should be called after GetWt + +SUBROUTINE GetBt(obj, elemsd, facetElemsd) + CLASS(TDGAlgorithm2_), INTENT(INOUT) :: obj + TYPE(ElemShapeData_), INTENT(IN) :: elemsd, facetElemsd + INTEGER(I4B) :: nrow, ncol + + nrow = obj%nrow + ncol = elemsd%nips + + obj%bt(1:nrow, 1:ncol) = MATMUL(obj%wmt(1:nrow, 1:nrow), & + elemsd%N(1:nrow, 1:ncol)) + + ! obj%bt(1:nrow, 1:ncol) = obj%bt(1:nrow, 1:ncol) * math%half + obj%bt(1:nrow, 1:ncol) = obj%bt(1:nrow, 1:ncol) + + obj%bt_right(1:nrow) = MATMUL(obj%wmt(1:nrow, 1:nrow), & + facetElemsd%N(1:nrow, 2)) + + ! obj%bt_right(1:nrow) = obj%bt_right(1:nrow) * math%half + obj%bt_right(1:nrow) = obj%bt_right(1:nrow) + +END SUBROUTINE GetBt + +!---------------------------------------------------------------------------- +! MakeZeros +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_MakeZeros +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_MakeZeros()" +#endif + +REAL(DFP), PARAMETER :: myzero = 0.0_DFP + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +obj%initialGuess_zero = obj%initialGuess.approxeq.myzero +obj%dis_zero = obj%dis.approxeq.myzero +obj%vel_zero = obj%vel.approxeq.myzero +obj%acc_zero = obj%acc.approxeq.myzero +obj%rhs_m_u1_zero = obj%rhs_m_u1.approxeq.myzero +obj%rhs_m_v1_zero = obj%rhs_m_v1.approxeq.myzero +obj%rhs_m_a1_zero = obj%rhs_m_a1.approxeq.myzero +obj%rhs_c_u1_zero = obj%rhs_c_u1.approxeq.myzero +obj%rhs_c_v1_zero = obj%rhs_c_v1.approxeq.myzero +obj%rhs_c_a1_zero = obj%rhs_c_a1.approxeq.myzero +obj%rhs_k_u1_zero = obj%rhs_k_u1.approxeq.myzero +obj%rhs_k_v1_zero = obj%rhs_k_v1.approxeq.myzero +obj%rhs_k_a1_zero = obj%rhs_k_a1.approxeq.myzero + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_MakeZeros + +!---------------------------------------------------------------------------- +! Deallocate +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Deallocate +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Deallocate()" +#endif + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +obj%isInit = .FALSE. + +obj%name = "TDG1" +obj%nrow = 0_I4B +obj%ncol = 0_I4B + +obj%initialGuess = 0.0_DFP +obj%initialGuess_zero = .TRUE. + +obj%dis = 0.0_DFP +obj%dis_zero = .TRUE. + +obj%vel = 0.0_DFP +obj%vel_zero = .TRUE. + +obj%acc = 0.0_DFP +obj%acc_zero = .TRUE. + +obj%mt = 0.0_DFP +obj%ct = 0.0_DFP +obj%bt = 0.0_DFP +obj%bt_right = 0.0_DFP +obj%wt = 0.0_DFP +obj%wmt = 0.0_DFP +obj%at = 0.0_DFP +obj%at_right = 0.0_DFP +obj%tat = 0.0_DFP +obj%kt = 0.0_DFP + +obj%rhs_m_u1 = 0.0_DFP +obj%rhs_m_v1 = 0.0_DFP +obj%rhs_m_a1 = 0.0_DFP +obj%rhs_c_u1 = 0.0_DFP +obj%rhs_c_v1 = 0.0_DFP +obj%rhs_c_a1 = 0.0_DFP +obj%rhs_k_u1 = 0.0_DFP +obj%rhs_k_v1 = 0.0_DFP +obj%rhs_k_a1 = 0.0_DFP + +obj%rhs_m_u1_zero = .TRUE. +obj%rhs_m_v1_zero = .TRUE. +obj%rhs_m_a1_zero = .TRUE. +obj%rhs_c_u1_zero = .TRUE. +obj%rhs_c_v1_zero = .TRUE. +obj%rhs_c_a1_zero = .TRUE. +obj%rhs_k_u1_zero = .TRUE. +obj%rhs_k_v1_zero = .TRUE. +obj%rhs_k_a1_zero = .TRUE. + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_Deallocate + +!---------------------------------------------------------------------------- +! Include errors +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE ConstructorMethods diff --git a/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@IOMethods.F90 b/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@IOMethods.F90 new file mode 100644 index 00000000..9e674e08 --- /dev/null +++ b/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@IOMethods.F90 @@ -0,0 +1,115 @@ +! 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 +! + +SUBMODULE(TDGAlgorithm2_Class) IOMethods +USE Display_Method, ONLY: Display +IMPLICIT NONE +CONTAINS + +!---------------------------------------------------------------------------- +! Display +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_Display +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_Display()" +#endif + +INTEGER(I4B) :: nrow, ncol + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +CALL Display(msg, unitno=unitno) +CALL Display(obj%isInit, "isInit: ", unitno=unitno) + +IF (.NOT. obj%isInit) THEN +#ifdef DEBUG_VER + CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif + RETURN +END IF + +CALL Display(obj%name, "name: ", unitno=unitno) +CALL Display(obj%nrow, "nrow: ", unitno=unitno) +CALL Display(obj%ncol, "ncol: ", unitno=unitno) + +nrow = obj%nrow +ncol = obj%ncol + +CALL Display(obj%initialGuess(1:nrow + 3), "initialGuess: ", & + unitno=unitno, advance="NO") +CALL Display(obj%initialGuess_zero(1:nrow + 3), "initialGuess_zero: ", & + unitno=unitno) + +CALL Display(obj%dis(1:nrow + 3), "dis: ", unitno=unitno, advance="NO") +CALL Display(obj%dis_zero(1:nrow + 3), "dis_zero: ", unitno=unitno) + +CALL Display(obj%vel(1:nrow + 3), "vel: ", unitno=unitno, advance="NO") +CALL Display(obj%vel_zero(1:nrow + 3), "vel_zero: ", unitno=unitno) + +CALL Display(obj%acc(1:nrow + 3), "acc: ", unitno=unitno, advance="NO") +CALL Display(obj%acc_zero(1:nrow + 3), "acc_zero: ", unitno=unitno) + +CALL Display(obj%mt(1:nrow, 1:ncol), "mt: ", unitno=unitno) +CALL Display(obj%ct(1:nrow, 1:ncol), "ct: ", unitno=unitno) +CALL Display(obj%kt(1:nrow, 1:ncol), "kt: ", unitno=unitno) + +CALL Display(obj%rhs_m_u1(1:nrow), "rhs_m_u1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_m_u1_zero(1:nrow), "rhs_m_u1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_m_v1(1:nrow), "rhs_m_v1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_m_v1_zero(1:nrow), "rhs_m_v1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_m_a1(1:nrow), "rhs_m_a1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_m_a1_zero(1:nrow), "rhs_m_a1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_c_u1(1:nrow), "rhs_c_u1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_c_u1_zero(1:nrow), "rhs_c_u1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_c_v1(1:nrow), "rhs_c_v1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_c_v1_zero(1:nrow), "rhs_c_v1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_c_a1(1:nrow), "rhs_c_a1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_c_a1_zero(1:nrow), "rhs_c_a1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_k_u1(1:nrow), "rhs_k_u1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_k_u1_zero(1:nrow), "rhs_k_u1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_k_v1(1:nrow), "rhs_k_v1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_k_v1_zero(1:nrow), "rhs_k_v1_zero: ", unitno=unitno) + +CALL Display(obj%rhs_k_a1(1:nrow), "rhs_k_a1: ", unitno=unitno, advance="NO") +CALL Display(obj%rhs_k_a1_zero(1:nrow), "rhs_k_a1_zero: ", unitno=unitno) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END] ') +#endif +END PROCEDURE obj_Display + +!---------------------------------------------------------------------------- +! Include errors +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE IOMethods diff --git a/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@TomlMethods.F90 b/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@TomlMethods.F90 new file mode 100644 index 00000000..857e5fab --- /dev/null +++ b/src/submodules/TDGAlgorithms/src/TDGAlgorithm2_Class@TomlMethods.F90 @@ -0,0 +1,98 @@ +! 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 +! + +SUBMODULE(TDGAlgorithm2_Class) TomlMethods +USE TomlUtility, ONLY: GetValue +USE tomlf, ONLY: toml_get => get_value + +IMPLICIT NONE + +CONTAINS + +!---------------------------------------------------------------------------- +! ImportFromToml +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ImportFromToml1 +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml1()" +#endif + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START] ') +#endif + +CALL obj%DEALLOCATE() + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END]') +#endif +END PROCEDURE obj_ImportFromToml1 + +!---------------------------------------------------------------------------- +! ImportFromToml +!---------------------------------------------------------------------------- + +MODULE PROCEDURE obj_ImportFromToml2 +#ifdef DEBUG_VER +CHARACTER(*), PARAMETER :: myName = "obj_ImportFromToml2()" +LOGICAL(LGT) :: isok +#endif + +TYPE(toml_table), ALLOCATABLE :: table +TYPE(toml_table), POINTER :: node +INTEGER(I4B) :: origin, stat + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[START]') +#endif + +CALL GetValue(table=table, afile=afile, filename=filename) + +node => NULL() +CALL toml_get(table, tomlName, node, origin=origin, requested=.FALSE., & + stat=stat) + +#ifdef DEBUG_VER +isok = ASSOCIATED(node) +CALL AssertError1(isok, myName, & + 'following error occured while reading the toml file :: & + &cannot find ['//tomlName//"] table in config.") +#endif + +CALL obj%ImportFromToml(table=node) + +node => NULL() +DEALLOCATE (table) + +#ifdef DEBUG_VER +CALL e%RaiseInformation(modName//'::'//myName//' - '// & + '[END]') +#endif +END PROCEDURE obj_ImportFromToml2 + +!---------------------------------------------------------------------------- +! Include errors +!---------------------------------------------------------------------------- + +#include "../../include/errors.F90" + +END SUBMODULE TomlMethods diff --git a/src/submodules/VectorField/src/VectorField_Class@BodySourceMethods.F90 b/src/submodules/VectorField/src/VectorField_Class@BodySourceMethods.F90 index 26ad2c6e..3c4bd8b9 100644 --- a/src/submodules/VectorField/src/VectorField_Class@BodySourceMethods.F90 +++ b/src/submodules/VectorField/src/VectorField_Class@BodySourceMethods.F90 @@ -110,8 +110,6 @@ args=args, n=elemsd%nsd) END DO - CALL Display(forceVecQuad(1:elemsd%nsd, 1:elemsd%nips), "forceVecQuad") - CALL FEVariable_Set( & obj=forceVar, val=forceVecQuad(1:elemsd%nsd, 1:elemsd%nips), & rank=TypeFEVariableVector, varType=TypeFEVariableSpace, scale=1.0_DFP, &