Gridap distributed

NEWS.md

@@ -8,17 +8,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+  - Extended support of `TriangulationPortion` to boundary and skeleton triangulations.  Since PR ...
   - Added `FESpaceWithLinearConstraints`. Since PR [#247](https://github.com/gridap/Gridap.jl/pull/247).
   - Added inner constructor to `CartesianDiscreteModel` allowing to build a model that represents a subgrid of
     a larger grid. Since PR [#245](https://github.com/gridap/Gridap.jl/pull/245).
 
 ### Changed
 
   - The part associated with the imposition of constraints in the `FESpace` interface has changed slightly. Since PR [#247](https://github.com/gridap/Gridap.jl/pull/247).
+  - Simplified the signature of `zero_free_values(::FESpace)`. Since PR ...
+  - Major refactoring in the `Assembler` interface. 
+    **Important change:** assembly-related functions take the data returned by functions like
+    `collect_cell_matrix` as it is. Example: the old user code `assemble_matrix(assembler,collect_cell_matrix(du,dv,terms)...)`
+    now is written simply as `assemble_matrix(assembler,collect_cell_matrix(du,dv,terms))`, i.e., the unpack of the last argument is not
+    used anymore.  With new new assembler interface, it is also possible to customize the assembly process
+    via a so-called `AssemblerStrategy` object. Since PR ...
   - Change the types of the sizes and partition fields of CartesianDescriptor to tuples instead of points. 
     Since PR [#246](https://github.com/gridap/Gridap.jl/pull/246).
 
-
 ## [0.9.2] - 2020-4-26
 
 ### Added

src/Algebra/Algebra.jl

@@ -13,6 +13,19 @@ using Test
 
 using Gridap.Helpers
 
+import Base: convert, size, getindex, show, count, *
+import LinearAlgebra: mul!
+import SparseArrays: nnz, nonzeros, nzrange, findnz, rowvals
+
+export allocate_vector
+export allocate_matrix
+export allocate_matrix_and_vector
+export allocate_in_domain
+export allocate_in_range
+export add_entries!
+export scale_entries!
+export muladd!
+
 export push_coo!
 export is_entry_stored
 export finalize_coo!
@@ -71,12 +84,18 @@ export get_vector
 export SparseMatrixCSR
 export SymSparseMatrixCSR
 
-include("MethodsAbstractMatrices.jl")
+include("AlgebraInterfaces.jl")
 
 include("SparseMatrices.jl")
 
 include("CompressedSparseMatrices.jl")
 
+include("SparseMatrixCSC.jl")
+
+include("SparseMatrixCSR.jl")
+
+include("SymSparseMatrixCSR.jl")
+
 include("NonlinearOperators.jl")
 
 include("NonlinearSolvers.jl")

src/Algebra/AlgebraInterfaces.jl

@@ -0,0 +1,135 @@
+
+function allocate_matrix end
+function allocate_matrix_and_vector end
+
+"""
+    allocate_vector(::Type{V},indices) where V
+
+Allocate a vector of type `V` indexable at the incides `indices`
+"""
+function allocate_vector(::Type{V},indices) where V
+  n = length(indices)
+  allocate_vector(V,n)
+end
+
+function allocate_vector(::Type{V},n::Integer) where V
+  T = eltype(V)
+  zeros(T,n)
+end
+
+"""
+    allocate_in_range(::Type{V},matrix) where V
+
+Allocate a vector of type `V` in the range of matrix `matrix`.
+"""
+function allocate_in_range(::Type{V},matrix) where V
+  n = size(matrix,1)
+  allocate_vector(V,n)
+end
+
+"""
+    allocate_in_domain(::Type{V},matrix) where V
+
+Allocate a vector of type `V` in the domain of matrix `matrix`.
+"""
+function allocate_in_domain(::Type{V},matrix) where V
+  n = size(matrix,2)
+  allocate_vector(V,n)
+end
+
+"""
+    fill_entries!(a,v)
+
+Fill the entries of array `a` with the value `v`. Returns `a`.
+"""
+function fill_entries!(a,v)
+  fill!(a,v)
+  a
+end
+
+"""
+    copy_entries!(a,b)
+
+Copy the entries of array `b` into array `a`. Returns `a`.
+"""
+function copy_entries!(a,b)
+  if a !== b
+    copyto!(a,b)
+  end
+  a
+end
+
+"""
+    add_entries!(a,b,combine=+)
+
+Perform the operation `combine` element-wise in the entries of arrays `a` and `b`
+and store the result in array `a`. Returns `a`.
+"""
+function add_entries!(a,b,combine=+)
+  @assert length(a) == length(b)
+  @inbounds for i in eachindex(a)
+    a[i] = combine(a[i],b[i])
+  end
+  a
+end
+
+"""
+    add_entry!(A,v,i,j,combine=+)
+
+Add an entry given its position and the operation to perform.
+"""
+function add_entry!(A,v,i::Integer,j::Integer,combine=+)
+  aij = A[i,j]
+  A[i,j] = combine(aij,v)
+end
+
+"""
+    add_entry!(A,v,i,combine=+)
+
+Add an entry given its position and the operation to perform.
+"""
+function add_entry!(A,v,i::Integer,combine=+)
+  ai = A[i]
+  A[i] = combine(ai,v)
+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)
+
+Matrix multiply a*b and add to result to c. Returns c.
+"""
+function muladd!(c,a,b)
+  _muladd!(c,a,b)
+  c
+end
+
+@static if VERSION >= v"1.3"
+  function _muladd!(c,a,b)
+    mul!(c,a,b,1,1)
+  end
+else
+  function _muladd!(c,a,b)
+    @assert length(c) == size(a,1)
+    @assert length(b) == size(a,2)
+    @inbounds for j in 1:size(a,2)
+      for i in 1:size(a,1)
+        c[i] += a[i,j]*b[j]
+      end
+    end
+  end
+end
+

src/Algebra/CompressedSparseMatrices.jl

@@ -1,34 +1,4 @@
 
-# CompressedSparseMatrices implementation contains:
-# 
-#   - Data types:
-#     + `SparseMatrixCSR`:Compressed Sparse Row (CSR) sparse matrix implementation with Bi-based indexing.
-#     + `SymSparseMatrixCSR`: Symmetric Compressed Sparse Row sparse matrix implementation with Bi-based indexing.
-# 
-#   - Procedures:
-#     + `sparsecsr`: Analogous method to [`sparse`](@ref) function to build a SparseMatrixCSR.
-#     + `symsparsecsr`: Analogous method to [`sparse`](@ref) function to build a SymSparseMatrixCSR.
-#     + `colvals`: Analogous method to [`rowvals`](@ref) to return `colvals` array.
-#     + `hasrowmajororder`: Return `true` if matrix values are ordered by row.
-#     + `hascolmajororder`: Return `true` if matrix values are ordered by column.
-#     + `getptr`: Return the pointer array of a SparseMatrix (`rowptr` or `colptr` depending on the SparseMatrix type)
-#     + `getindices`: Return the indices array of a SparseMatrix (`rowval` or `colval` depending on the SparseMatrix type)
-# 
-#   - Overloaded procedures:
-#     + `*`: SparseMatrix-Vector product.
-#     + `mul!`: SparseMatrix-Vector product.
-#     + `nnz`: Return the number of stored (filled) elements in a sparse array.
-#     + `nonzeros`: Return `nzval` array.
-#     + `nzrange`: Return the range of indices for a particular row or column (Depending on the SparseMatrix type)
-#     + `findnz`: Return a tuple (I, J, V) where I and J are the row and column indices.
-#     + `rowvals`: Return row indices or raises an error (Depending on the SparseMatrix type)
-#     + `convert`: Type conversion
-#     + `count`: Count the number of elements in nonzeros(S) for which predicate pred returns true. 
-
-import Base: convert, size, getindex, show, count, *
-import LinearAlgebra: mul!
-import SparseArrays: nnz, nonzeros, nzrange, findnz, rowvals
-
 """
     hasrowmajororder(::Type{AbstractSparseMatrix})
 
@@ -157,7 +127,3 @@ function fill_entries!(A::AbstractSparseMatrix{Tv,Ti},v::Number) where {Tv,Ti}
   nonzeros(A) .= convert(Tv,v)
 end
 
-
-include("SparseMatrixCSC.jl")
-include("SparseMatrixCSR.jl")
-include("SymSparseMatrixCSR.jl")

src/Algebra/LinearSolvers.jl

@@ -20,7 +20,7 @@ get_vector(op::AffineOperator) = op.vector
 
 function residual!(b::AbstractVector,op::AffineOperator,x::AbstractVector)
   mul!(b,op.matrix,x)
-  b .-=  op.vector
+  add_entries!(b,op.vector,-)
   b
 end
 
@@ -33,13 +33,14 @@ function jacobian(op::AffineOperator,x::AbstractVector)
   op.matrix
 end
 
-function zero_initial_guess(::Type{T},op::AffineOperator) where T
-  n = size(op.matrix,2)
-  zeros(T,n)
+function zero_initial_guess(op::AffineOperator)
+  x = allocate_in_domain(typeof(op.vector),op.matrix)
+  fill_entries!(x,zero(eltype(x)))
+  x
 end
 
 function allocate_residual(op::AffineOperator,x::AbstractVector)
-  similar(x,eltype(x),length(op.vector))
+  allocate_in_range(typeof(op.vector),op.matrix)
 end
 
 function allocate_jacobian(op::AffineOperator,x::AbstractVector)
@@ -151,12 +152,12 @@ function solve!(x::AbstractVector,
                 ls::LinearSolver,
                 op::NonlinearOperator,
                 cache::Nothing)
-  x .= zero(eltype(x))
+  fill_entries!(x,zero(eltype(x)))
   b = residual(op, x)
   A = jacobian(op, x)
   ss = symbolic_setup(ls, A)
   ns = numerical_setup(ss,A)
-  broadcast!(*,b,b,-1)
+  scale_entries!(b,-1)
   solve!(x,ns,b)
   LinearSolverCache(A,b,ns)
 end
@@ -165,13 +166,13 @@ function solve!(x::AbstractVector,
                 ls::LinearSolver,
                 op::NonlinearOperator,
                 cache)
-  x .= zero(eltype(x))
+  fill_entries!(x,zero(eltype(x)))
   b = cache.b
   A = cache.A
   ns = cache.ns
   residual!(b, op, x)
   numerical_setup!(ns,A)
-  broadcast!(*,b,b,-1)
+  scale_entries!(b,-1)
   solve!(x,ns,b)
   cache
 end
@@ -238,7 +239,8 @@ end
 function solve!(
   x::AbstractVector,ns::LUNumericalSetup,b::AbstractVector)
   y = ns.factors\b # the allocation of y can be avoided
-  x .= y
+  copy_entries!(x,y)
+  x
 end
 
 """
@@ -265,5 +267,6 @@ end
 
 function solve!(
   x::AbstractVector,ns::BackslashNumericalSetup,b::AbstractVector)
-   copyto!(x, ns.A\b)
+  copy_entries!(x, ns.A\b)
 end
+

src/Algebra/NLSolvers.jl

@@ -54,9 +54,9 @@ function _solve_nr!(x,A,b,dx,ns,nls,op)
   for nliter in 1:nls.max_nliters
 
     # Solve linearized problem
-    broadcast!(*,b,b,-1)
+    scale_entries!(b,-1)
     solve!(dx,ns,b)
-    broadcast!(+,x,x,dx)
+    add_entries!(x,dx)
 
     # Check convergence for the current residual
     residual!(b, op, x)
@@ -157,7 +157,7 @@ function _nlsolve_with_updated_cache!(x,nls,op,cache)
   end
   r = nlsolve(df,x;linsolve=linsolve!,kwargs...)
   cache.result = r
-  x[:] .= r.zero
+  copy_entries!(x,r.zero)
 end
 
 function _new_nlsolve_cache(x0,nls,op)