Merge branch 'philippe'

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.2"
 manifest_format = "2.0"
-project_hash = "4ff3722b2aaa4099ffdfd4d945b58147c681fc55"
+project_hash = "0de19c1262248bf2550740d15ddd11277591182e"
 
 [[deps.AbstractFFTs]]
 deps = ["ChainRulesCore", "LinearAlgebra"]
@@ -61,7 +61,7 @@ uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
 version = "0.4.2"
 
 [[deps.Cairo_jll]]
-deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
+deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
 git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
 version = "1.16.1+1"
@@ -325,9 +325,9 @@ version = "0.21.0+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "fb83fbe02fe57f2c068013aa94bcdf6760d3a7a7"
+git-tree-sha1 = "d3b3624125c1474292d0d8ed0f65554ac37ddb23"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.74.0+1"
+version = "2.74.0+2"
 
 [[deps.Graphics]]
 deps = ["Colors", "LinearAlgebra", "NaNMath"]
@@ -545,9 +545,9 @@ version = "1.42.0+0"
 
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
+git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.16.1+1"
+version = "1.16.1+2"
 
 [[deps.Libmount_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -580,6 +580,12 @@ git-tree-sha1 = "2ce8695e1e699b68702c03402672a69f54b8aca9"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
 version = "2022.2.0+0"
 
+[[deps.MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.10"
+
 [[deps.Makie]]
 deps = ["Animations", "Base64", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Distributions", "DocStringExtensions", "FFMPEG", "FileIO", "FixedPointNumbers", "Formatting", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MakieCore", "Markdown", "Match", "MathTeXEngine", "MiniQhull", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "Printf", "Random", "RelocatableFolders", "Serialization", "Showoff", "SignedDistanceFields", "SnoopPrecompile", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
 git-tree-sha1 = "0462fdf13cfad9fa6f23a2fb41da38e7a2aa4908"

Project.toml

@@ -7,6 +7,7 @@ version = "0.2.1"
 EspyInsideFunction = "d6b55740-5cb4-11ec-1cb9-93db9984c461"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/make.jl

@@ -1,8 +1,7 @@
-using Documenter, EspyInsideFunction
+using Documenter, Muscade
 
 push!(LOAD_PATH,"../src")
 makedocs(sitename ="Muscade.jl",
-#        modules   = [EspyInsideFunction],
         format    = Documenter.HTML(prettyurls = false),
         pages     = ["index.md"],
         source  = "src",

docs/src/index.md

@@ -1,32 +1,25 @@
+# Muscade.jl
+
 ```@meta
 CurrentModule = Muscade
 ```
 
-# Muscade.jl
+## [Purpose](@id purpose)
 
-Documentation for [Muscade.jl](https://github.com/SINTEF/Muscade.jl).
+[Muscade.jl](https://github.com/SINTEF/Muscade.jl) is a framework to support the rapid development of domain-specific, multiphysics optimization-FEM models.
 
-```@index
-```
+Muscade is a framework: it does not provide finite element types modeling any specific domain of physics: rather, it facilitates their development.  Using automatic differentiation and automatic result extraction, the code needed to implement new elements is minimized.  The API to create new elements is quite classic (given degrees of freedom, compute their duals, aka. residuals): Muscade is not a *modeling language* like for example UFL which provides a formalism in which a "well posed problem" can be defined with domain, differential equations and boundary conditions.
 
-## [Some chapter](@id some-chapter)
+Degrees of freedom (dofs) in Muscade have a `field`: application developers can introduce new fields like (`:tx` for translation in `x` direction, `c` for solute concentration, etc...), and associate a scale to them, to aid solvers handle unknowns of differents orders of magnitudes.  This facilitates the implementation of multiphysical models.
 
-The following is an example of annotated code
+Muscade handles optimization-FEM problems, that is, optimization problems constrained by equilibrium of the FEM model.
 
-```jldoctest EspyDemo; output = false
-1+1==2
-# output
-true
-```
+## [Reference manual](@id reference)
 
-```julia
-some code
+```@index
 ```
 
-
-## Reference manual
-
-
 ```@autodocs
 Modules = [Muscade]
+Order   = [:module,:constant,:type,:function,:macro, :type]
 ```

src/core/Adiff.jl

@@ -4,11 +4,11 @@ using   Printf
 
 
 ## Type and construction
-AA = AbstractArray
-AV = AbstractVector
-SV = SVector  
-SA = SArray 
-SM = SMatrix
+const AA = AbstractArray
+const AV = AbstractVector
+const SV = SVector  
+const SA = SArray 
+const SM = SMatrix
 # Types
 # P precedence 
 # N number of partials 
@@ -99,11 +99,12 @@ value{P}(a::∂ℝ{P,N,R}) where{P,N,R   } = a.x
 value{P}(a::R        ) where{P  ,R<:ℝ} = a
 value{P}(a::AA{R}    ) where{P  ,R   } = value{P}.(a)
 
-# no ∂{P}(a) syntax: in case a does not contain adiffs 
+# ∂{P}(a) is handled as ∂{P,1}(a) and returns a scalar 
 ∂{P,N}(a::     ∂ℝ{P,N,R} ) where{  P,N,R   } = a.dx
 ∂{P,N}(a::            R  ) where{  P,N,R<:ℝ} = SV{  N,R}(zero(R)    for i=1:N      )
 ∂{P,N}(a::SV{M,∂ℝ{P,N,R}}) where{M,P,N,R   } = SM{M,N,R}(a[i].dx[j] for i=1:M,j∈1:N) # ∂(a,x)[i,j] = ∂a[i]/∂x[j]
 ∂{P,N}(a::SV{M,       R }) where{M,P,N,R   } = SM{M,N,R}(zero(R)    for i=1:M,j=1:N)
+∂{P  }(a::            R  ) where{  P,  R<:ℝ} = zero(R)
 ∂{P  }(a::     ∂ℝ{P,1,R} ) where{  P,  R   } = a.dx[1]
 ∂{P  }(a::SV{N,∂ℝ{P,1,R}}) where{  P,N,R   } = SV{  N,R}(a[i].dx[1] for i=1:N     ) # ∂(a,x)[i]    = ∂a[i]/∂x
 #∂{P,N}(a::SA{M,∂ℝ{P,N,R}}) where{M,P,N,R}  = SA{(M...,N),R}(a[i].dx[j] for i∈eachindex(a),j∈1:N) # ∂(a,x)[i,...,j] = ∂a[i,...]/∂x[j]

src/core/Assemble.jl

@@ -195,13 +195,13 @@ function asmvec!(asm,dofgr,dis)
     Λ,X,U,A  = indexedstate(dofgr)      # create a state of indices into the group - with zeros for modeldofs not in group
     for (ieletyp,di) ∈ enumerate(dis.dis)
         nΛ,nX,nU,nA = gradientstructure(dofgr,di) # number of dofs of each class in the gradient returned by an element
-        iΛ,iX,iU,iA = gradientpartition(nΛ,nX,nU,nA)  # indices into said gradient
+        iΛ,iX,iU,iA = gradientpartition(nΛ,nX,nU,nA)  # indices into said gradient TODO type unstable, barrier function
         asm[ieletyp] = zeros(𝕫,nΛ+nX+nU+nA,length(di.index)) # asm[ieletyp][idof,iele] (its a view)
         for (iele,index) ∈ enumerate(di.index)
-            asm[ieletyp][iΛ,iele] = nonzeros(Λ[index.X])  
-            asm[ieletyp][iX,iele] = nonzeros(X[index.X])
-            asm[ieletyp][iU,iele] = nonzeros(U[index.U])
-            asm[ieletyp][iA,iele] = nonzeros(A[index.A])
+            asm[ieletyp][iΛ,iele] .= nonzeros(Λ[index.X])  
+            asm[ieletyp][iX,iele] .= nonzeros(X[index.X])
+            asm[ieletyp][iU,iele] .= nonzeros(U[index.U])
+            asm[ieletyp][iA,iele] .= nonzeros(A[index.A])
         end
     end
     return 𝕣1(undef,getndof(dofgr))
@@ -319,30 +319,45 @@ function assemblesequential!(out,asm,dis,eleobj,state::State{Nxder,Nuder},γ,dbg
     end
 end
 
-#### addin and zero!
+#### addtoarray and zero!
 function zero!(out::DenseArray)
     out .= 0
 end
 function zero!(out::AbstractSparseArray)
     out.nzval .= 0
 end
 
-function addin!(out::DenseArray,asm,iele,a) 
-    for (i,ai) ∈ enumerate(a)
-        j = asm[i,iele]
-        if j≠0
-            out[j]+=ai
+function add_value!(out::𝕣1,asm,iele,a::SVector{M,∂ℝ{P,N,𝕣}},ias) where{P,N,M}
+    for (iasm,ia) ∈ enumerate(ias)
+        iout = asm[iasm,iele]
+        if iout≠0
+            out[iout]+=a[ia].x
         end
     end
 end   
-function addin!(out::AbstractSparseArray,asm,iele,a)
-    for (i,ai) ∈ enumerate(a)
-        j = asm[i,iele]
-        if j≠0
-            out.nzval[j]+=ai
+function add_value!(out::𝕣1,asm,iele,a::SVector{M,𝕣},ias) where{M}
+    for (iasm,ia) ∈ enumerate(ias)
+        iout = asm[iasm,iele]
+        if iout≠0
+            out[iout]+=a[ia]
         end
     end
-end
+end   
+add_value!(out,asm,iele,a) = add_value!(out,asm,iele,a,eachindex(a)) 
+struct add_∂!{P} end 
+function add_∂!{P}(out::Array,asm,iele,a::SVector{M,∂ℝ{P,N,R}},i1as,i2as) where{P,N,R,M}
+    for (i1asm,i1a) ∈ enumerate(i1as), (i2asm,i2a) ∈ enumerate(i2as)
+        iasm = i1asm+length(i1as)*(i2asm-1)
+        iout = asm[iasm,iele]
+        if iout≠0
+            out[iout]+=a[i1a].dx[i2a]  
+        end
+    end
+end  
+add_∂!{P}(out::SparseMatrixCSC,args...) where{P}                      = add_∂!{P}(out.nzval,args...)
+add_∂!{P}(out::Array,asm,iele,a::SVector{M,R},args...) where{P,M,R}   = nothing
+add_∂!{P}(out::Array,asm,iele,a::SVector{M,∂ℝ{P,N,R}}) where{P,N,R,M} = add_∂!{P}(out,asm,iele,a,1:M,1:N)
+
 
 ###### scaled functions