Merge pull request #171 from OpenSEMBA/dev

Dev

Author: lmdiazangulo

doc/smbjson.md

@@ -576,9 +576,9 @@ Performs a loop integral along on the contour of the surface reference in the `e
 
 + If it is a point or a volume, `direction` must be present.
 + For an oriented line, `direction` is optional and the orientation of the line will be used as default.
-+ For an oriented surface, `direction` is also optional, and if not value is given it is assumed to be the normal of the surface.
++ For an oriented surface, `direction` is also optional, and if not value is given it is assumed to be the the positive axis of the surface normal.
 
-Due to Ampere's law, the loop integral of the magnetic field is equal to the total electric current passing through the surfaces. `[field]`, can be `electric` or `magnetic`. Defaults to `electric`, which gives the total current passing through the surface.
+Due to Ampere's law, the loop integral of the magnetic field is equal to the total electric current passing through the surfaces. `[field]`, can be `electric` or `magnetic`. Defaults to `electric`, which returns the total current passing through the surface.
 
 In the following example `elementId` points an element describing a single oriented surface, therefore `direction` does not need to be stated explicitly.
 
@@ -704,7 +704,7 @@ An example of a planewave propagating towards $\hat{z}$ and polarized in the $+\
 
 This object represents a time-varying vector field applied along an oriented line with the same orientation of the line. Therefore, the `elementIds` within must contain only elements of type `cell` with `intervals` describing a collection of oriented lines. Additionally, it may contain:
 
-+ `[field]` with a `electric` or `magnetic` label which indicates the vector field which will be applied. If not present, it defaults to `electric`.
++ `[field]` with a `current` label which indicates the vector field which will be applied. If not present, it defaults to `current`.
 + `[hardness]` with `soft` or `hard` label. A `soft` hardness indicated that the magnitude will be **added** to the field this situation is typical for a waveport. `hard` sources mean that the field is **substituted** by the value established by the `magnitudeFile`, which for an electric field `nodalSource` would be equivalent to a `pec` material if the magnitude is zero.
 
 **Example:**

src_json_parser/nfdetypes_extension.F90

@@ -673,9 +673,7 @@ elemental logical function curr_field_src_eq(a, b)
          (a%n_C2P == b%n_C2P) .and. &
          (a%nombre == b%nombre) .and. &
          (a%isElec .eqv. b%isElec) .and. &
-         (a%isMagnet  .eqv. b%isMagnet) .and. &
-         (a%isCurrent .eqv. b%isCurrent) .and. &
-         (a%isField .eqv. b%isField) .and. &
+         (a%isHard .eqv. b%isHard) .and. &
          (a%isInitialValue .eqv. b%isInitialValue)
    end function curr_field_src_eq
 

src_json_parser/smbjson.F90

@@ -832,30 +832,36 @@ function readField(jns) result(res)
          type(Curr_Field_Src) :: res
          type(json_value), pointer :: jns, entry
          integer, dimension(:), allocatable :: elementIds
+         character (len=BUFSIZE) :: nodalSourceName
          type(coords_scaled), dimension(:), allocatable :: coordsFromLinels
 
-         select case (this%getStrAt(jns, J_FIELD))
-          case (J_FIELD_ELECTRIC)
-            res%isField = .true.
+         select case (this%getStrAt(jns, J_FIELD, default=J_FIELD_CURRENT))
+          case (J_FIELD_CURRENT)
             res%isElec = .true.
-            res%isMagnet = .false.
-            res%isCurrent = .false.
-          case (J_FIELD_MAGNETIC)
-            res%isField = .true.
-            res%isElec = .false.
-            res%isMagnet = .true.
-            res%isCurrent = .false.
           case default
             write(error_unit, *) 'Error reading current field source. Field label not recognized.'
          end select
+         
+         select case (this%getStrAt(jns, J_SRC_NS_HARDNESS, default=J_SRC_NS_HARDNESS_SOFT))
+          case (J_SRC_NS_HARDNESS_SOFT)
+            res%isHard = .false.
+          case (J_SRC_NS_HARDNESS_HARD)
+            res%isHard = .true.
+          case default
+            write(error_unit, *) 'Error reading current field source. Hardness label not recognized.'
+          end select
+         
          res%isInitialValue = .false.
+
          res%nombre = trim(adjustl(this%getStrAt(jns, J_SRC_MAGNITUDE_FILE)))
 
+         nodalSourceName = this%getStrAt(jns, J_NAME, default=' ')
+
          allocate(res%c1P(0))
          res%n_C1P = 0
 
          elementIds = this%getIntsAt(jns, J_ELEMENTIDS)
-         call cellRegionsToScaledCoords(res%C2P,  this%mesh%getCellRegions(elementIds) )
+         call cellRegionsToScaledCoords(res%C2P,  this%mesh%getCellRegions(elementIds), nodalSourceName)
          res%n_C2P = size(res%C2p)
 
       end function
@@ -1225,7 +1231,7 @@ function readBlockProbe(bp) result(res)
          res%j2  = cs(1)%ye
          res%k1  = cs(1)%zi
          res%k2  = cs(1)%ze
-         res%nml = cs(1)%Or
+         res%nml = abs(cs(1)%Or)
 
          res%outputrequest = trim(adjustl(this%getStrAt(bp, J_NAME)))
          call setDomain(res, this%getDomain(bp, J_PR_DOMAIN))

src_json_parser/smbjson_labels.F90

@@ -164,6 +164,11 @@ module smbjson_labels_mod
    character (len=*), parameter :: J_SRC_PW_THETA = "theta"
    character (len=*), parameter :: J_SRC_PW_PHI = "phi"
 
+   ! type(NodalSource)
+   character (len=*), parameter :: J_SRC_NS_HARDNESS = "hardness"
+   character (len=*), parameter :: J_SRC_NS_HARDNESS_SOFT = "soft"
+   character (len=*), parameter :: J_SRC_NS_HARDNESS_HARD = "hard"
+
    ! --- probe types
    character (len=*), parameter :: J_PROBES = "probes"
 

src_main_pub/nfde_types.F90

@@ -678,8 +678,7 @@ MODULE NFDETypes
       CHARACTER (LEN=BUFSIZE) :: nombre
       INTEGER (KIND=4) :: n_C1P = 0
       INTEGER (KIND=4) :: n_C2P = 0
-      LOGICAL :: isElec, isMagnet
-      LOGICAL :: isCurrent, isField, isInitialValue
+      LOGICAL :: isElec, isHard, isInitialValue
    END TYPE Curr_Field_Src
    !------------------------------------------------------------------------------
    ! Definicin de las Nodal Source global

src_main_pub/preprocess_geom.F90

@@ -2896,8 +2896,8 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
             !
             sgg%NodalSource(conta1)%fichero%name = trim (adjustl(this%nodsrc%NodalSource(i)%nombre))
             sgg%NodalSource(conta1)%isElec = this%nodsrc%NodalSource(i)%isElec
-            sgg%NodalSource(conta1)%IsHard = this%nodsrc%NodalSource(i)%isField
-            sgg%NodalSource(conta1)%IsInitialValue = this%nodsrc%NodalSource(i)%IsInitialValue
+            sgg%NodalSource(conta1)%IsHard = this%nodsrc%NodalSource(i)%isHard
+            sgg%NodalSource(conta1)%IsInitialValue = this%nodsrc%NodalSource(i)%IsInitialValue 
          endif
          !
          tama2 = this%nodsrc%NodalSource(i)%n_c1P
@@ -2929,7 +2929,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         !bug OLD 181214 sl_4_20mm_gli.nfde. Fuente nodal electrica embebida en pec y nodal magnetica en pmc se ignoraran sean hard or soft
                         MEDIO = sggmiEx (i1, j1, k1)
                         valido=.true.
-                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         ! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !ELSE
@@ -2944,7 +2944,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         ! COMENTADO 250816 PQ DA UN ERROR JUISTO CUANDO CAE LA FUENTE EN UN CORTE MPI. HABRIA QU TOCA LA CASUISTICA DE punto_s%ZE = Min (this%nodsrc%NodalSource(i)%c1P(ii)%ZE, Max(BoundingBox%ZI, BoundingBox%ZE  )) PERO NO LO HE QUERIDO HACER
                         !!!MEDIO = sggmiHx (i1, j1, k1)
                         !!!valido=.true.
-                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !!!!  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         !!!! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !!!!ELSE
@@ -2961,7 +2961,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                      IF (punto_s%yc /= 0) THEN
                         MEDIO = sggmiEy (i1, j1, k1)
                         valido=.true.
-                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         ! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !ELSE
@@ -2976,7 +2976,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         !
                         !!!MEDIO = sggmiHy (i1, j1, k1)
                         !!!valido=.true.
-                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !!!!  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         !!!! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !!!!ELSE
@@ -2993,7 +2993,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                      IF (punto_s%zc /= 0) THEN
                         MEDIO = sggmiEz (i1, j1, k1)
                         valido=.true.
-                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         ! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !ELSE
@@ -3008,7 +3008,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         !
                         !!!MEDIO = sggmiHz (i1, j1, k1)
                         !!!valido=.true.
-                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !!!!  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         !!!! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !!!!ELSE
@@ -3103,7 +3103,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                      IF (punto_s%xc /= 0) THEN
                         MEDIO = sggmiEx (i1, j1, k1)
                         valido=.true.
-                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         ! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !ELSE
@@ -3119,7 +3119,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         !
                         !!!MEDIO = sggmiHx (i1, j1, k1)
                         !!!valido=.true.
-                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !!!!  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         !!!! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !!!!ELSE
@@ -3136,7 +3136,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                      IF (punto_s%yc /= 0) THEN
                         MEDIO = sggmiEy (i1, j1, k1)
                         valido=.true.
-                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         ! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !ELSE
@@ -3151,7 +3151,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         !
                         !!!MEDIO = sggmiHy (i1, j1, k1)
                         !!!valido=.true.
-                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !!!!  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         !!!! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !!!!ELSE
@@ -3168,7 +3168,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                      IF (punto_s%zc /= 0) THEN
                         MEDIO = sggmiEz (i1, j1, k1)
                         valido=.true.
-                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         ! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !ELSE
@@ -3183,7 +3183,7 @@ SUBROUTINE read_geomData (sgg,sggMtag,tag_numbers, sggMiNo,sggMiEx,sggMiEy,sggMi
                         !
                         !!!MEDIO = sggmiHz (i1, j1, k1)
                         !!!valido=.true.
-                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isField) THEN
+                        !!!!IF ( .NOT. this%nodsrc%NodalSource(i)%isHard) THEN
                         !!!!  VALIDO = (sgg%Med(MEDIO)%Is%Dielectric) .OR. (sgg%Med(MEDIO)%Is%EDispersive) .OR. &
                         !!!! & (sgg%Med(MEDIO)%Is%MDispersive)
                         !!!!ELSE

src_pyWrapper/pyWrapper.py

@@ -10,6 +10,7 @@
 from vtk.util.numpy_support import vtk_to_numpy
 from itertools import product
 import copy
+import matplotlib.pyplot as plt
 
 DEFAULT_SEMBA_FDTD_PATH = '/build/bin/semba-fdtd'
 
@@ -194,6 +195,19 @@ def _positionStrToTwoCells(pos_str: str):
     def _getFieldAndDirection(tag: str):
         return tag[1], tag[2]
 
+class ExcitationFile():
+    def __init__(self, excitation_filename):
+        if isinstance(excitation_filename, os.PathLike):
+            self.filename = excitation_filename.as_posix()
+        else:
+            self.filename = excitation_filename
+        assert os.path.isfile(self.filename)
+
+        self.data = pd.read_csv(self.filename, sep='\\s+', names=['time', 'value'])
+
+    def __getitem__(self, key):
+        return self.data[key]
+    
 
 class FDTD():
     def __init__(self, input_filename, path_to_exe=None,
@@ -335,6 +349,19 @@ def getSolvedProbeFilenames(self, probe_name):
 
         return sorted(probeFiles)
 
+    def getExcitationFile(self, excitation_file_name):
+        file_extensions =('*.1.exc',)
+        excitationFile = []
+        for ext in file_extensions:
+            newExcitationFile = [x for x in glob.glob(ext) if re.match(excitation_file_name, x)]
+            excitationFile.extend(newExcitationFile)
+        
+        if ((len(excitationFile)) != 1):
+            raise "Unexpected number of excitation Files found: {}".format(excitationFile)
+
+        return excitationFile
+
+
     def getVTKMap(self):
         current_path = os.getcwd()
         folders = [item for item in os.listdir(

test/pyWrapper/test_full_system.py

@@ -450,6 +450,37 @@ def testCanExecuteFDTDFromFolderWithSpacesAndCanProcessAdditionalArguments(tmp_p
     solver.run()
     assert (Probe(solver.getSolvedProbeFilenames("outside")[0]) is not None)
     assert (solver.getVTKMap()[0] is not None)
+    
+
+def test_nodal_source(tmp_path):
+    fn = CASES_FOLDER + "nodalSource/nodalSource.fdtd.json"
+    assert (os.path.isfile(fn))
+    solver = FDTD(fn, path_to_exe=SEMBA_EXE, run_in_folder=tmp_path)
+    solver.run()
+    
+    resistanceBulkProbe = Probe( \
+    solver.getSolvedProbeFilenames("Bulk probe Resistance")[0])
+    nodalBulkProbe = Probe( \
+        solver.getSolvedProbeFilenames("Bulk probe Nodal Source")[0])
+    excitation = ExcitationFile( \
+        excitation_filename=solver.getExcitationFile("predefinedExcitation")[0])
+
+    # For debugging.
+    # plt.figure()
+    # plt.plot(resistanceBulkProbe['time'].to_numpy(), 
+    #         resistanceBulkProbe['current'].to_numpy(), label='BP Current@resistance')
+    # plt.plot(excitation.data['time'].to_numpy(), 
+    #         excitation.data['value'].to_numpy(), label='excited current')
+    # plt.plot(nodalBulkProbe['time'].to_numpy(),
+    #         -nodalBulkProbe['current'].to_numpy(), label='BP Current@nodal source')
+    # plt.legend()
+
+    exc = np.interp(nodalBulkProbe['time'].to_numpy(), 
+                    excitation.data['time'].to_numpy(), 
+                    excitation.data['value'].to_numpy())
+    assert np.corrcoef(exc, -nodalBulkProbe['current'])[0,1] > 0.999
+    assert np.corrcoef(-nodalBulkProbe['current'], resistanceBulkProbe['current'])[0,1] > 0.998
+
 
 def testCanAssignSameSurfaceImpedanceToMultipleGeometries(tmp_path):
     fn = CASES_FOLDER + 'multipleAssigments/multipleSurfaceImpedance.fdtd.json'

test/pyWrapper/test_pyWrapper.py

@@ -212,4 +212,4 @@ def test_casemaker_sphere_rcs_case(tmp_path):
         case_name + ".fdtd.json",
         path_to_exe=SEMBA_EXE)
     solver.run()
-    assert solver.hasFinishedSuccessfully()
+    assert solver.hasFinishedSuccessfully()

test/smbjson/test_read_airplane.F90

@@ -65,9 +65,7 @@ function expectedProblemDescription() result (expected)
       allocate(expected%nodSrc%NodalSource(1))
       expected%nodSrc%NodalSource(1)%nombre = "gauss.exc"
       expected%nodSrc%NodalSource(1)%isElec = .true.
-      expected%nodSrc%NodalSource(1)%isMagnet = .false.
-      expected%nodSrc%NodalSource(1)%isCurrent = .false.
-      expected%nodSrc%NodalSource(1)%isField = .true.
+      expected%nodSrc%NodalSource(1)%isHard = .false.
       expected%nodSrc%NodalSource(1)%isInitialValue = .false.
       allocate(expected%nodSrc%NodalSource(1)%c1P(0))
       allocate(expected%nodSrc%NodalSource(1)%c2P(1))
@@ -79,7 +77,7 @@ function expectedProblemDescription() result (expected)
       expected%nodSrc%NodalSource(1)%c2P(1)%Ye = 30
       expected%nodSrc%NodalSource(1)%c2P(1)%Zi = 39
       expected%nodSrc%NodalSource(1)%c2P(1)%Ze = 46
-      expected%nodSrc%NodalSource(1)%c2P(1)%tag = ''
+      expected%nodSrc%NodalSource(1)%c2P(1)%tag = 'nodalSource'
       expected%nodSrc%NodalSource(1)%c2P(1)%xc = 0.0
       expected%nodSrc%NodalSource(1)%c2P(1)%yc = 0.0
       expected%nodSrc%NodalSource(1)%c2P(1)%zc = 1.0