Merge pull request #1050 from Antoinemarteau/master

Fix AttachDirichletMap type stability

Author: JordiManyer

NEWS.md

@@ -5,6 +5,10 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+### Fixed
+
+- Fixed #974, an error when weak form is real but unknown vector is complex. Since PR[#1050](https://github.com/gridap/Gridap.jl/pull/1050).
+
 ## [0.18.7] - 2024-10-8
 
 ### Added

src/Arrays/AlgebraMaps.jl

@@ -35,6 +35,7 @@ function evaluate!(cache,k::MulAddMap,a,b,c)
   setsize!(cache,size(c))
   d = cache.array
   copyto!(d,c)
+  iszero(k.α) && isone(k.β) && return d
   mul!(d,a,b,k.α,k.β)
   d
 end

src/CellData/AttachDirichlet.jl

@@ -70,13 +70,14 @@ function Arrays.return_cache(k::AttachDirichletMap,matvec::Tuple,vals,mask)
 end
 
 function Arrays.evaluate!(cache,k::AttachDirichletMap,matvec::Tuple,vals,mask)
+  mat, vec = matvec
   if mask
-    mat, vec = matvec
     vec_with_bcs = evaluate!(cache,k.muladd,mat,vals,vec)
-    (mat, vec_with_bcs)
   else
-    matvec
+    identity_muladd = MulAddMap(0,1)
+    vec_with_bcs = evaluate!(cache,identity_muladd ,mat,vals,vec)
   end
+  (mat, vec_with_bcs)
 end
 
 function Arrays.return_value(k::AttachDirichletMap,mat,vals,mask)

src/Fields/ArrayBlocks.jl

@@ -1281,6 +1281,7 @@ function evaluate!(cache,k::MulAddMap,a::ArrayBlock,b::ArrayBlock,c::ArrayBlock)
   _setsize_mul!(c1,a,b)
   d = evaluate!(c2,unwrap_cached_array,c1)
   copyto!(d,c)
+  iszero(k.α) && isone(k.β) && return d
   mul!(d,a,b,k.α,k.β)
   d
 end
@@ -1517,4 +1518,4 @@ function Base.show(io::IO,o::ArrayBlockView)
 end
 
 LinearAlgebra.diag(a::MatrixBlockView) = view(a.array.array, diag(a.block_map))
-LinearAlgebra.diag(a::MatrixBlock) = view(a.array,diag(CartesianIndices(a.array)))
+LinearAlgebra.diag(a::MatrixBlock) = view(a.array,diag(CartesianIndices(a.array)))

test/GridapTests/issue_974.jl

@@ -0,0 +1,31 @@
+using Gridap
+
+T = ComplexF64
+
+# Operator creating a linear system where the matrix and RHS are Float64 but the unknown vector is ComplexF64.
+domain = (0, 1, 0, 1)
+partition = (4, 4)
+model = CartesianDiscreteModel(domain, partition)
+
+order = 1
+reffe = ReferenceFE(lagrangian, Float64, order)
+
+V1 = TestFESpace(model, reffe; conformity=:H1, vector_type=Vector{T})
+V2 = TestFESpace(model, reffe; conformity=:H1, vector_type=Vector{T})
+
+U1 = TrialFESpace(V1)
+U2 = TrialFESpace(V2)
+
+Y = MultiFieldFESpace([V1, V2])
+X = MultiFieldFESpace([U1, U2])
+
+degree = 2 * order
+Ω = Triangulation(model)
+dΩ = Measure(Ω, degree)
+
+# Project constant 1 into both spaces.
+a((u1, u2), (v1, v2)) = ∫(v1 * u1)dΩ + ∫(v2 * u2)dΩ
+l((v1, v2)) = ∫(v1 + v2)dΩ
+
+op = AffineFEOperator(a, l, X, Y)
+uh1, uh2 = solve(op)