Merge pull request #74 from WIAS-PDELib/feature/coupling

Use sparse vector in periodic coupling

Author: pjaap

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableFEM"
 uuid = "a722555e-65e0-4074-a036-ca7ce79a4aed"
-version = "1.3.0"
+version = "1.3.1"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]
@@ -29,7 +29,7 @@ DiffResults = "1"
 DocStringExtensions = "0.8,0.9"
 ExampleJuggler = "2.2.1"
 ExplicitImports = "1"
-ExtendableFEMBase = "1"
+ExtendableFEMBase = "1.3.0"
 ExtendableGrids = "1.10.3"
 ExtendableSparse = "1.5.3"
 ForwardDiff = "0.10.35,1"

src/ExtendableFEM.jl

@@ -67,7 +67,7 @@ using LinearSolve: LinearSolve, LinearProblem, UMFPACKFactorization, deleteat!,
     init, solve
 using Printf: Printf, @printf, @sprintf
 using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, findnz, nnz,
-    nzrange, rowvals, sparse
+    nzrange, rowvals, sparse, SparseVector
 using SparseDiffTools: SparseDiffTools, ForwardColorJacCache,
     forwarddiff_color_jacobian!, matrix_colors
 using Symbolics: Symbolics

src/helper_functions.jl

@@ -244,8 +244,9 @@ function _get_periodic_coupling_matrix(
     # to be sure
     fill!(fe_vector.entries, 0.0)
 
-    # FE vector for interpolation
-    fe_vector_target = FEVector(FES)
+    # FE target vector for interpolation with sparse entries
+    fe_vector_target = FEVector(FES; entries = SparseVector{Float64, Int64}(FES.ndofs, Int64[], Float64[]))
+
 
     # resulting sparse matrix
     n = length(fe_vector.entries)
@@ -358,6 +359,13 @@ function _get_periodic_coupling_matrix(
         append!(searchareas, view(faces_to, 1:nfaces_to))
     end
 
+    # throw error if no search area had been found for a bface
+    for source in 1:num_sources(searchareas)
+        if num_targets(searchareas, source) == 0
+            throw("bface $source has no valid search area on the opposite side of the grid. Are from/to regions and give_opposite! function correct?")
+        end
+    end
+
     # loop over boundary face indices: we need this index for dofs_on_boundary
     for i_boundary_face in 1:n_boundary_faces
 
@@ -372,7 +380,8 @@ function _get_periodic_coupling_matrix(
                 end
 
                 # reset
-                fill!(fe_vector_target.entries, 0.0)
+                empty!(fe_vector_target.entries.nzind)
+                empty!(fe_vector_target.entries.nzval)
 
                 # activate one entry
                 fe_vector.entries[local_dof] = 1.0
@@ -394,7 +403,7 @@ function _get_periodic_coupling_matrix(
                 fe_vector.entries[local_dof] = 0.0
 
                 # set entries
-                for (i, target_entry) in enumerate(fe_vector_target.entries)
+                for (i, target_entry) in zip(findnz(fe_vector_target.entries)...)
                     if abs(target_entry) > sparsity_tol
                         result[local_dof, i] = target_entry
                     end