Merge pull request #688 from gridap/release_0.17

Release 0.17

NEWS.md

@@ -4,6 +4,23 @@ 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).
 
+## [0.17.0] - 2021-10-22
+
+### Added
+
+- Aliases `Interior`, `Boundary`, `Skeleton`, and `Interface` for the `Triangulation`, `BoundaryTriangulation`, `SkeletonTriangulation`, and `InterfaceTriangulation` constructors. Since PR [#662](https://github.com/gridap/Gridap.jl/pull/662).
+- Function `create_pvtk_file` for exporting results in `pvtu` format. Since PR [#685](https://github.com/gridap/Gridap.jl/pull/685).
+
+### Changed
+
+- Major refactoring in the `Triangulation` interface to properly support the solution of PDEs defined on domains of different dimension. The major change from the user perspective is that `Triangulation` objects can be used both to integrate the weak form (as before) but also to define FE spaces (except for unfitted triangulations obviously). It is still possible to define FE spaces from `DiscreteModels`, but it is safer and more idiomatic (closer to the math notation) to use `Triangulation` objects from now on. Since PR [#662](https://github.com/gridap/Gridap.jl/pull/662).
+- Changes in assembly interface to allow optimization when assembling matrices and vectors simultaneously. Since PR [#685](https://github.com/gridap/Gridap.jl/pull/685).
+
+### Removed
+
+- `BoundaryDiscreteModel`, `RestrictedDiscreteMdeol`, `RestrictedTriangulation`, `TriangulationStyle`, `BackgroundTriangulation`, `SubTriangulation`, `get_cell_to_bgcell`, `get_cell_ref_map`, `get_background_triangulation`, and `have_compatible_domains`. Since PR [#662](https://github.com/gridap/Gridap.jl/pull/662). 
+- Functions `scale_entries!` and `fill_entries!`. Replaced by Julia functions `LinearAlgebra.rmul!` and `LinearAlgebra.fillstored!`. Since PR [#680](https://github.com/gridap/Gridap.jl/pull/680).
+
 ## [0.16.5] - 2021-09-08
 
 ### Added

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <[email protected]>", "Francesc Verdugo <[email protected]>"]
-version = "0.16.5"
+version = "0.17.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -43,9 +43,9 @@ LineSearches = "7.0.1"
 NLsolve = "4.3.0"
 NearestNeighbors = "0.4.8"
 QuadGK = "2.3.1, 2.4"
-SparseMatricesCSR = "0.6"
 StaticArrays = "0.12.1, 1.0"
-WriteVTK = "1.7, 1.8"
+WriteVTK = "=1.11.0"
+SparseMatricesCSR = "0.6.4"
 julia = "1.3"
 
 [extras]

src/Algebra/Algebra.jl

@@ -26,7 +26,6 @@ export allocate_matrix_and_vector
 export allocate_in_domain
 export allocate_in_range
 export add_entries!
-export scale_entries!
 export muladd!
 export nz_counter
 export nz_allocation
@@ -45,7 +44,6 @@ export is_entry_stored
 export finalize_coo!
 export sparse_from_coo
 export add_entry!
-export fill_entries!
 export copy_entries!
 export allocate_coo_vectors
 export push_coo!

src/Algebra/AlgebraInterfaces.jl

@@ -61,17 +61,6 @@ function allocate_in_domain(::Type{V},matrix) where V
   allocate_vector(V,n)
 end
 
-"""
-    fill_entries!(a,v)
-
-Fill the entries of array `a` with the value `v`. Returns `a`.
-For sparse matrices it only fills the non-zero entries.
-"""
-function fill_entries!(a,v)
-  fill!(a,v)
-  a
-end
-
 """
     copy_entries!(a,b)
 
@@ -196,19 +185,6 @@ end
   A
 end
 
-"""
-    scale_entries!(a,v)
-
-Scale the entries of array `a` with the value `v`. Returns `a`.
-"""
-function scale_entries!(a,b)
-  @inbounds for i in eachindex(a)
-    a[i] = b*a[i]
-  end
-  a
-end
-
-# Base.mul!
 
 """
     muladd!(c,a,b)
@@ -280,7 +256,7 @@ end
 # We can also do a loop and update
 # the entries of c
 #
-#    fill_entries!(c,0)
+#    fill!(c,0) or LinearAlgebra.fillstored!(c,0)
 #    add_entry!(c,v,i,j)
 #    add_entries!(c,vs,is,js)
 #
@@ -290,6 +266,20 @@ struct DoNotLoop end
 LoopStyle(::Type) = DoNotLoop()
 LoopStyle(::T) where T = LoopStyle(T)
 
+# By default process, matrix and vector separately
+# but, in some situations, create_from_nz of the vector
+# can reuse data from the one computed in
+# create_from_nz for the matrix (e.g., GridapDistributed)
+function create_from_nz(a,b)
+  A = create_from_nz(a)
+  B = create_from_nz(b)
+  A,B
+end
+# See comment above for create_from_nz. The same applies here
+# for nz_allocation.
+function nz_allocation(a,b)
+  nz_allocation(a),nz_allocation(b)
+end
 
 # For dense arrays
 
@@ -342,34 +332,34 @@ SparseMatrixBuilder(a::SparseMatrixBuilder) = a
 
 get_array_type(::SparseMatrixBuilder{T}) where T = T
 
-mutable struct SparseMatrixCounter{T,A}
+mutable struct CounterCOO{T,A}
   nnz::Int
   axes::A
-  function SparseMatrixCounter{T}(axes::A) where {T,A<:NTuple{2,AbstractUnitRange}}
+  function CounterCOO{T}(axes::A) where {T,A<:NTuple{2,AbstractUnitRange}}
     nnz = 0
     new{T,A}(nnz,axes)
   end
 end
 
-LoopStyle(::Type{<:SparseMatrixCounter}) = Loop()
+LoopStyle(::Type{<:CounterCOO}) = Loop()
 
-@inline function add_entry!(::Function,a::SparseMatrixCounter{T},v,i,j) where T
+@inline function add_entry!(::Function,a::CounterCOO{T},v,i,j) where T
   if is_entry_stored(T,i,j)
     a.nnz = a.nnz + 1
   end
   a
 end
 
-struct CooAllocation{T,A,B,C}
-  counter::SparseMatrixCounter{T,A}
+struct AllocationCOO{T,A,B,C}
+  counter::CounterCOO{T,A}
   I::B
   J::B
   V::C
 end
 
-LoopStyle(::Type{<:CooAllocation}) = Loop()
+LoopStyle(::Type{<:AllocationCOO}) = Loop()
 
-@inline function add_entry!(::typeof(+),a::CooAllocation{T},::Nothing,i,j) where T
+@inline function add_entry!(::typeof(+),a::AllocationCOO{T},::Nothing,i,j) where T
   if is_entry_stored(T,i,j)
     a.counter.nnz = a.counter.nnz + 1
     k = a.counter.nnz
@@ -379,7 +369,7 @@ LoopStyle(::Type{<:CooAllocation}) = Loop()
   nothing
 end
 
-@inline function add_entry!(::typeof(+),a::CooAllocation{T},v,i,j) where T
+@inline function add_entry!(::typeof(+),a::AllocationCOO{T},v,i,j) where T
   if is_entry_stored(T,i,j)
     a.counter.nnz = a.counter.nnz + 1
     k = a.counter.nnz
@@ -391,20 +381,20 @@ end
 end
 
 #function nz_counter(::Type{T},axes) where T<:AbstractSparseMatrix
-#  SparseMatrixCounter{T}(axes)
+#  CounterCOO{T}(axes)
 #end
 
 function nz_counter(::SparseMatrixBuilder{T},axes) where T<:AbstractSparseMatrix
-  SparseMatrixCounter{T}(axes)
+  CounterCOO{T}(axes)
 end
 
-function nz_allocation(a::SparseMatrixCounter{T}) where T
-  counter = SparseMatrixCounter{T}(a.axes)
+function nz_allocation(a::CounterCOO{T}) where T
+  counter = CounterCOO{T}(a.axes)
   I,J,V = allocate_coo_vectors(T,a.nnz)
-  CooAllocation(counter,I,J,V)
+  AllocationCOO(counter,I,J,V)
 end
 
-function create_from_nz(a::CooAllocation{T}) where T
+function create_from_nz(a::AllocationCOO{T}) where T
   m,n = map(length,a.counter.axes)
   finalize_coo!(T,a.I,a.J,a.V,m,n)
   sparse_from_coo(T,a.I,a.J,a.V,m,n)
@@ -466,13 +456,7 @@ function copy_entries!(a::T,b::T) where T<:AbstractSparseMatrix
   end
 end
 
-function fill_entries!(A::AbstractSparseMatrix,v)
-  nonzeros(A) .= v
-  A
-end
-
 function allocate_coo_vectors(
    ::Type{<:AbstractSparseMatrix{Tv,Ti}},n::Integer) where {Tv,Ti}
   (zeros(Ti,n), zeros(Ti,n), zeros(Tv,n))
 end
-

src/Algebra/LinearSolvers.jl

@@ -43,7 +43,7 @@ end
 
 function zero_initial_guess(op::AffineOperator)
   x = allocate_in_domain(typeof(op.vector),op.matrix)
-  fill_entries!(x,zero(eltype(x)))
+  fill!(x,zero(eltype(x)))
   x
 end
 
@@ -174,12 +174,12 @@ function solve!(x::AbstractVector,
                 ls::LinearSolver,
                 op::NonlinearOperator,
                 cache::Nothing)
-  fill_entries!(x,zero(eltype(x)))
+  fill!(x,zero(eltype(x)))
   b = residual(op, x)
   A = jacobian(op, x)
   ss = symbolic_setup(ls, A)
   ns = numerical_setup(ss,A)
-  scale_entries!(b,-1)
+  rmul!(b,-1)
   solve!(x,ns,b)
   LinearSolverCache(A,b,ns)
 end
@@ -188,13 +188,13 @@ function solve!(x::AbstractVector,
                 ls::LinearSolver,
                 op::NonlinearOperator,
                 cache)
-  fill_entries!(x,zero(eltype(x)))
+  fill!(x,zero(eltype(x)))
   b = cache.b
   A = cache.A
   ns = cache.ns
   residual!(b, op, x)
   numerical_setup!(ns,A)
-  scale_entries!(b,-1)
+  ruml!(b,-1)
   solve!(x,ns,b)
   cache
 end

src/Algebra/NLSolvers.jl

@@ -54,7 +54,7 @@ function _solve_nr!(x,A,b,dx,ns,nls,op)
   for nliter in 1:nls.max_nliters
 
     # Solve linearized problem
-    scale_entries!(b,-1)
+    rmul!(b,-1)
     solve!(dx,ns,b)
     x .+= dx
 
@@ -183,5 +183,3 @@ function _update_nlsolve_cache!(cache,x0,op)
   numerical_setup!(ns,j0)
   NLSolversCache(f0,j0,df,ns,nothing)
 end
-
-

src/Arrays/AlgebraMaps.jl

@@ -60,8 +60,3 @@ end
 function evaluate!(cache,k::TouchEntriesMap,A,v,i)
   add_entries!(+,A,nothing,i)
 end
-
-
-
-
-

src/Arrays/AppendedArrays.jl

@@ -111,6 +111,19 @@ function lazy_map(::typeof(evaluate),::Type{T},b::Fill,a::AppendedArray...) wher
   end
 end
 
+function lazy_map(::typeof(evaluate),::Type{T},b::Fill,a::AppendedArray) where T
+  f = b.value
+  ra = lazy_map(f,a.a)
+  rb = lazy_map(f,a.b)
+  lazy_append(ra,rb)
+end
+
+function lazy_map(k::Reindex{<:LazyArray},j_to_i::AppendedArray)
+  a = lazy_map(k,j_to_i.a)
+  b = lazy_map(k,j_to_i.b)
+  lazy_append(a,b)
+end
+
 # In that case we don't implement it for ::Type{T} variant since we want to
 # avoid to run type inference since the first entry in a can have non-concrete
 # eltype

src/Arrays/Arrays.jl

@@ -25,9 +25,6 @@ import Base: getindex, setindex!
 import Base: similar
 import Base: IndexStyle
 
-import Gridap.Algebra: scale_entries!
-import Gridap.Algebra: fill_entries!
-
 # CachedArray
 
 export CachedArray

src/Arrays/PosNegReindex.jl

@@ -211,3 +211,4 @@ all_neg(a::PosNegPartition) = length(a.ipos_to_i)==0
 #   i = j_to_i[j]
 #   i_to_v[i]=v
 # end
+

src/CellData/AttachDirichlet.jl

@@ -83,11 +83,9 @@ end
   cm, cv = cache
   if mask
     vec_with_bcs = evaluate!(cm,*,mat,vals)
-    scale_entries!(vec_with_bcs,-1)
+    rmul!(vec_with_bcs,-1)
   else
     vec_with_bcs = evaluate!(cv,ZeroVectorMap(),mat)
   end
   (mat, vec_with_bcs)
 end
-
-

src/CellData/CellData.jl

@@ -32,7 +32,7 @@ import Gridap.Geometry: num_cells
 import Gridap.Geometry: get_triangulation
 
 import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part
-import LinearAlgebra: det, tr, cross, dot, ⋅
+import LinearAlgebra: det, tr, cross, dot, ⋅, rmul!
 import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj
 
 export gradient, ∇