Merge branch 'master' of github.com:gridap/Gridap.jl into moment-base…

…d-reffes

.github/workflows/ci.yml

@@ -42,7 +42,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           file: lcov.info
           verbose: true

NEWS.md

@@ -5,9 +5,21 @@ 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.18.8] - 2024-12-2
+
+### Added
+
+- Added get_dof_value_type for FESpacesWithLinearConstraints. Since PR[#1062](https://github.com/gridap/Gridap.jl/pull/1062).
+- Added Xiao-Gimbutas quadratures for simplices. Since PR[#1058](https://github.com/gridap/Gridap.jl/pull/1058).
+- Small improvements of the documentation of `Gridap.TensorValues`. Since PR[#1051](https://github.com/gridap/Gridap.jl/pull/1051).
+
 ### 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).
+- Fixed issue where barycentric refinement rule in 3D would not produce oriented meshes. Since PR[#1055](https://github.com/gridap/Gridap.jl/pull/1055).
+
+### Changed
+- Optimized MonomialBasis low-level functions. Since PR[#1059](https://github.com/gridap/Gridap.jl/pull/1059).
 
 ## [0.18.7] - 2024-10-8
 

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 authors = ["Santiago Badia <[email protected]>", "Francesc Verdugo <[email protected]>", "Alberto F. Martin <[email protected]>"]
-version = "0.18.7"
+version = "0.18.8"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

benchmark/Project.toml

@@ -2,3 +2,4 @@
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

benchmark/benchmarks.jl

@@ -12,3 +12,4 @@ end
 const SUITE = BenchmarkGroup()
 
 @include_bm SUITE "bm_assembly"
+@include_bm SUITE "bm_monomial_basis"

benchmark/bm/bm_monomial_basis.jl

@@ -0,0 +1,190 @@
+module bm_monomial_basis
+
+using PkgBenchmark, BenchmarkTools
+using Gridap
+using Gridap.Polynomials
+using Gridap.TensorValues
+using StaticArrays
+
+################################################
+# src/Polynomials/MonomialBasis.jl: _set_value_!
+################################################
+
+gradient_type = Gridap.Fields.gradient_type
+
+_set_value! = Gridap.Polynomials._set_value!
+
+function set_value_driver(f,T,D,x,n)
+  k = 1
+  s = one(T)
+  for i in 1:n
+    k = f(x,s,k)
+  end
+end
+
+function set_value_benchmarkable(D, T, V, n)
+  C = num_indep_components(V)
+  x = zeros(V,n*C)
+  return @benchmarkable set_value_driver($_set_value!,$T,$D,$x,$n)
+end
+
+##################################################
+# src/Polynomials/ModalC0Bases.jl: _set_value_mc0!
+##################################################
+
+_set_value_mc0! = Gridap.Polynomials._set_value_mc0!
+
+function set_value_mc0_driver(f,T,D,x,n)
+  k = 1
+  s = one(T)
+  for i in 1:n
+    k = f(x,s,k,2)
+  end
+end
+
+function set_value_mc0_benchmarkable(D, T, V, n)
+  C = num_indep_components(V)
+  x = zeros(V,2*n*C)
+  return @benchmarkable set_value_mc0_driver($_set_value_mc0!,$T,$D,$x,$n)
+end
+
+###################################################
+# src/Polynomials/MonomialBasis.jl: _set_gradient!
+###################################################
+
+ _set_gradient! = Gridap.Polynomials. _set_gradient!
+
+function set_gradient_driver(f,T,D,V,x,n)
+  k = 1
+  s = VectorValue{D,T}(ntuple(_->one(T),D))
+  for i in 1:n
+    k = f(x,s,k,V)
+  end
+end
+
+function set_gradient_benchmarkable(D, T, V, n)
+  C = num_indep_components(V)
+  G = gradient_type(V, zero(Point{D,T}))
+  x = zeros(G,n*C);
+  return @benchmarkable set_gradient_driver($_set_gradient!,$T,$D,$V,$x,$n)
+end
+
+#####################################################
+# src/Polynomials/ModalC0Bases.jl: _set_gradient_mc0!
+#####################################################
+
+ _set_gradient_mc0! = Gridap.Polynomials. _set_gradient_mc0!
+
+function set_gradient_mc0_driver(f,T,D,V,x,n)
+  k = 1
+  s = VectorValue{D,T}(ntuple(_->one(T),D))
+  for i in 1:n
+    k = f(x,s,k,1,V)
+  end
+end
+
+function set_gradient_mc0_benchmarkable(D, T, V, n)
+  C = num_indep_components(V)
+  G = gradient_type(V, zero(Point{D,T}))
+  x = zeros(G,n*C);
+  return @benchmarkable set_gradient_mc0_driver($_set_gradient_mc0!,$T,$D,$V,$x,$n)
+end
+
+#################################################
+# src/Polynomials/MonomialBasis.jl: _evaluate_1d!
+#################################################
+
+_evaluate_1d! = Gridap.Polynomials._evaluate_1d!
+
+function evaluate_1d_driver(f,order,D,v,x_vec)
+  for x in x_vec
+    f(v,x,order,D)
+  end
+end
+
+function evaluate_1d_benchmarkable(D, T, V, n)
+  n = Integer(n/50)
+  order = num_indep_components(V)
+  v = zeros(D,order+1);
+  x = rand(MVector{n,T})
+  return @benchmarkable evaluate_1d_driver($_evaluate_1d!,$order,$D,$v,$x)
+end
+
+################################################
+# src/Polynomials/MonomialBasis.jl:_gradient_1d!
+################################################
+
+_gradient_1d! = Gridap.Polynomials._gradient_1d!
+
+function gradient_1d_driver(f,order,D,v,x_vec)
+  for x in x_vec
+    f(v,x,order,D)
+  end
+end
+
+function gradient_1d_benchmarkable(D, T, V, n)
+  n = Integer(n/10)
+  order = num_indep_components(V)
+  v = zeros(D,order+1);
+  x = rand(MVector{n,T})
+  return @benchmarkable gradient_1d_driver($_gradient_1d!,$order,$D,$v,$x)
+end
+
+################################################
+# src/Polynomials/MonomialBasis.jl:_hessian_1d!
+################################################
+
+_hessian_1d! = Gridap.Polynomials._hessian_1d!
+
+function hessian_1d_driver(f,order,D,v,x_vec)
+  for x in x_vec
+    f(v,x,order,D)
+  end
+end
+
+function hessian_1d_benchmarkable(D, T, V, n)
+  n = Integer(n/10)
+  order = num_indep_components(V)
+  v = zeros(D,order+1);
+  x = rand(MVector{n,T})
+  return @benchmarkable hessian_1d_driver($_hessian_1d!,$order,$D,$v,$x)
+end
+
+#####################
+# benchmarkable suite
+#####################
+
+const SUITE = BenchmarkGroup()
+
+const benchmarkables = (
+  set_value_benchmarkable,
+  set_value_mc0_benchmarkable,
+  set_gradient_benchmarkable,
+  set_gradient_mc0_benchmarkable,
+  evaluate_1d_benchmarkable,
+  gradient_1d_benchmarkable,
+  hessian_1d_benchmarkable
+)
+
+const dims=(1, 2, 3, 5, 8)
+const n = 3000
+const T = Float64
+
+for benchable in benchmarkables
+  for D in dims
+    TV = [
+      VectorValue{D,T},
+      TensorValue{D,D,T,D*D},
+      SymTensorValue{D,T,Integer(D*(D+1)/2)},
+      SymTracelessTensorValue{D,T,Integer(D*(D+1)/2)}
+    ]
+
+    for V in TV
+      if V == SymTracelessTensorValue{1,T,1} continue end # no dofs
+      name = "monomial_basis_$(D)D_$(V)_$(benchable)"
+      SUITE[name] = benchable(D, T, V, n)
+    end
+  end
+end
+
+end # module

docs/src/TensorValues.md

@@ -1,9 +1,151 @@
+# Gridap.TensorValues
+
 ```@meta
 CurrentModule = Gridap.TensorValues
 ```
 
-# Gridap.TensorValues
+This module provides the abstract interface `MultiValue` representing tensors
+that are also `Number`s, along with concrete implementations for the following
+tensors:
+- 1st order [`VectorValue`](@ref),
+- 2nd order [`TensorValue`](@ref),
+- 2nd order and symmetric [`SymTensorValue`](@ref),
+- 2nd order, symmetric and traceless [`SymTracelessTensorValue`](@ref),
+- 3rd order [`ThirdOrderTensorValue`](@ref),
+- 4th order and symmetric [`SymFourthOrderTensorValue`](@ref).
+
+## Generalities
+
+The main feature of this module is that the provided types do not extend from `AbstractArray`, but from `Number`!
+
+This allows one to work with them as if they were scalar values in broadcasted operations on arrays of `VectorValue` objects (also for `TensorValue` or `MultiValue` objects). For instance, one can perform the following manipulations:
+```julia
+# Assign a VectorValue to all the entries of an Array of VectorValues
+A = zeros(VectorValue{2,Int}, (4,5))
+v = VectorValue(12,31)
+A .= v # This is possible since  VectorValue <: Number
 
-```@autodocs
-Modules = [TensorValues,]
+# Broadcasting of tensor operations in arrays of TensorValues
+t = TensorValue(13,41,53,17) # creates a 2x2 TensorValue
+g = TensorValue(32,41,3,14) # creates another 2x2 TensorValue
+B = fill(t,(1,5))
+C = inner.(g,B) # inner product of g against all TensorValues in the array B
+@show C
+# C = [2494 2494 2494 2494 2494]
 ```
+
+To create a [`::MultiValue`](@ref) tensor from components, these should be given
+as separate arguments or all gathered in a `tuple`. The order of the arguments
+is the order of the linearized Cartesian indices of the corresponding array
+(order of the `Base.LinearIndices` indices):
+```julia
+using StaticArrays
+t = TensorValue( (1, 2, 3, 4) )
+ts= convert(SMatrix{2,2,Int}, t)
+@show ts
+# 2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
+#  1  3
+#  2  4
+t2[1,2] == t[1,2] == 3 # true
+```
+For symmetric tensor types, only the independent components should be given, see
+[`SymTensorValue`](@ref), [`SymTracelessTensorValue`](@ref) and [`SymFourthOrderTensorValue`](@ref).
+
+A `MultiValue` can be created from an `AbstractArray` of the same size. If the
+`MultiValue` type has internal constraints (e.g. symmetries), ONLY the required
+components are picked from the array WITHOUT CHECKING if the given array
+did respect the constraints:
+```julia
+SymTensorValue( [1 2; 3 4] )          # -> SymTensorValue{2, Int64, 3}(1, 2, 4)
+SymTensorValue( SMatrix{2}(1,2,3,4) ) # -> SymTensorValue{2, Int64, 3}(1, 3, 4)
+```
+
+`MultiValue`s can be converted to static and mutable arrays types from
+`StaticArrays.jl` using `convert` and [`mutable`](@ref), respectively.
+
+## Tensor types
+
+The following concrete tensor types are currently implemented:
+
+```@docs
+VectorValue
+TensorValue
+SymTensorValue
+SymTracelessTensorValue
+ThirdOrderTensorValue
+SymFourthOrderTensorValue
+```
+
+### Abstract tensor types
+
+```@docs
+MultiValue
+AbstractSymTensorValue
+```
+
+## Interface
+
+The tensor types implement methods for the following `Base` functions: `getindex`, `length`, `size`, `rand`, `zero`, `real`, `imag` and `conj`.
+
+`one` is also implemented in particular cases: it is defined for second
+and fourth order tensors. For second order, it returns the identity tensor `δij`,
+except `SymTracelessTensorValue` that does not implement `one`. For fourth order symmetric tensors, see [`one`](@ref).
+
+Additionally, the tensor types expose the following interface:
+
+```@docs
+num_components
+mutable
+Mutable
+num_indep_components
+indep_comp_getindex
+indep_components_names
+change_eltype
+
+inner
+dot
+double_contraction
+outer
+```
+
+### Other type specific interfaces
+
+#### For square second order tensors
+
+```@docs
+det
+inv
+symmetric_part
+```
+
+#### For first order tensors
+
+```@docs
+diagonal_tensor
+```
+
+#### For second and third order tensors
+
+```@docs
+tr
+```
+
+#### For first and second order tensors
+
+```@docs
+norm
+meas
+```
+
+#### For `VectorValue` of length 2 and 3
+
+```@docs
+cross
+```
+
+#### For second order non-traceless and symmetric fourth order tensors
+
+```@docs
+one
+```
+