Merge branch 'master' of github.com:Jai-Tushar/devGridap.jl into hho

Project.toml

@@ -24,6 +24,7 @@ PolynomialBases = "c74db56a-226d-5e98-8bb0-a6049094aeea"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

src/CellData/CellQuadratures.jl

@@ -1,22 +1,25 @@
 
 # Quadrature rule
+const QET=Union{<:Quadrature,<:VectorBlock{<:Quadrature}}
+const PET=Union{<:AbstractArray{<:Point},<:VectorBlock{<:AbstractArray{<:Point}}}
+const WET=Union{<:AbstractArray{<:Real},<:VectorBlock{<:AbstractArray{<:Real}}}
 
 """
 """
 struct CellQuadrature{DDS,IDS} <: CellDatum
-  cell_quad::AbstractArray{<:Quadrature}
-  cell_point::AbstractArray{<:AbstractArray{<:Point}}
-  cell_weight::AbstractArray{<:AbstractArray{<:Real}}
+  cell_quad::AbstractArray{<:QET}
+  cell_point::AbstractArray{<:PET}
+  cell_weight::AbstractArray{<:WET}
   trian::Triangulation
   data_domain_style::DDS
   integration_domain_style::IDS
 end
 
 # Old constructor (for backward compatibility)
 function CellQuadrature(
-  cell_quad::AbstractArray{<:Quadrature},
-  cell_point::AbstractArray{<:AbstractArray{<:Point}},
-  cell_weight::AbstractArray{<:AbstractArray{<:Real}},
+  cell_quad::AbstractArray{<:QET},
+  cell_point::AbstractArray{<:PET},
+  cell_weight::AbstractArray{<:WET},
   trian::Triangulation,
   data_domain_style::DomainStyle)
   CellQuadrature(cell_quad,cell_point,cell_weight,trian,data_domain_style,PhysicalDomain())

src/FESpaces/Assemblers.jl

@@ -474,22 +474,15 @@ function collect_cell_matrix_and_vector(
   biform::DomainContribution,liform::DomainContribution)
 
   matvec, mat, vec = _pair_contribution_when_possible(biform,liform)
-  matvecdata = _collect_cell_matvec(trial,test,matvec)
-  matdata = collect_cell_matrix(trial,test,mat)
-  vecdata = collect_cell_vector(test,vec)
-  (matvecdata, matdata, vecdata)
+  _collect_cell_matrix_and_vector(trial, test, matvec, mat, vec)
 end
 
 function collect_cell_matrix_and_vector(
   trial::FESpace,test::FESpace,
   biform::DomainContribution,liform::DomainContribution,uhd::FEFunction)
 
   matvec, mat, vec = _pair_contribution_when_possible(biform,liform,uhd)
-
-  matvecdata = _collect_cell_matvec(trial,test,matvec)
-  matdata = collect_cell_matrix(trial,test,mat)
-  vecdata = collect_cell_vector(test,vec)
-  (matvecdata, matdata, vecdata)
+  _collect_cell_matrix_and_vector(trial, test, matvec, mat, vec)
 end
 
 function _pair_contribution_when_possible(biform,liform)
@@ -513,21 +506,35 @@ end
 
 function _pair_contribution_when_possible(biform,liform,uhd)
   _matvec, _mat, _vec = _pair_contribution_when_possible(biform,liform)
-  matvec = DomainContribution()
-  mat = DomainContribution()
-  for (trian,t) in _matvec.dict
+  matvec, mat = _attach_dirichlet(_matvec, _mat, uhd)
+  matvec, mat, _vec
+end
+
+function _collect_cell_matrix_and_vector(trial, test, matvec, mat, vec)
+  matvecdata = _collect_cell_matvec(trial,test,matvec)
+  matdata = collect_cell_matrix(trial,test,mat)
+  vecdata = collect_cell_vector(test,vec)
+  (matvecdata, matdata, vecdata)
+end
+
+function _attach_dirichlet(matvec, mat, uhd)
+  _matvec = DomainContribution()
+  _mat = DomainContribution()
+  for (trian,t) in matvec.dict
     cellvals = get_cell_dof_values(uhd,trian)
     cellmask = get_cell_is_dirichlet(uhd,trian)
-    matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
+    _matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
   end
-  for (trian,t) in _mat.dict
+  for (trian,t) in mat.dict
     cellvals = get_cell_dof_values(uhd,trian)
     cellmask = get_cell_is_dirichlet(uhd,trian)
-    matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
+    _matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
   end
-  matvec, mat, _vec
+  _matvec, _mat
 end
 
+
+
 # allow linear forms like `l(v) = 0
 
 function collect_cell_vector(test::FESpace,l::Number)

src/Fields/ArrayBlocks.jl

@@ -21,7 +21,7 @@
 # uh [c][lf][np]
 
 
-struct ArrayBlock{A,N}
+struct ArrayBlock{A,N} 
   array::Array{A,N}
   touched::Array{Bool,N}
   function ArrayBlock(array::Array{A,N},touched::Array{Bool,N}) where {A,N}

src/Geometry/AppendedTriangulations.jl

@@ -129,11 +129,11 @@ function get_glue(t::AppendedTriangulation,::Val{D}) where D
   lazy_append(a,b)
 end
 
-function lazy_append(a::FaceToFaceGlue,b::FaceToFaceGlue)
+function lazy_append(a::FaceToFaceGlue{Dc},b::FaceToFaceGlue{Dc}) where Dc
   tface_to_mface = lazy_append(a.tface_to_mface,b.tface_to_mface)
   tface_to_mface_map = lazy_append(a.tface_to_mface_map,b.tface_to_mface_map)
   mface_to_tface = nothing
-  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
+  FaceToFaceGlue(Dc,tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
 function get_cell_shapefuns(trian::AppendedTriangulation)

src/Geometry/BoundaryTriangulations.jl

@@ -204,7 +204,7 @@ function get_glue(trian::BoundaryTriangulation,::Val{D},::Val{D}) where D
   face_s_q = lazy_map(linear_combination,face_to_q_vertex_coords,face_to_shapefuns)
   tface_to_mface_map = face_s_q
   mface_to_tface = nothing
-  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
+  FaceToFaceGlue(num_cell_dims(trian),tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
 function get_facet_normal(trian::BoundaryTriangulation, boundary_trian_glue::FaceToCellGlue)
@@ -252,9 +252,13 @@ function push_normal(invJt,n)
   end
 end
 
-function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation,glue)
+function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation)
   d = num_cell_dims(trian)
-  cell_grid = get_grid(get_background_model(trian.trian))
+  _compute_face_to_q_vertex_coords_body(d,get_background_model(trian.trian),trian.glue)
+end
+
+function _compute_face_to_q_vertex_coords_body(d,model,glue)
+  cell_grid = get_grid(model)
   polytopes = map(get_polytope, get_reffes(cell_grid))
   cell_to_ctype = glue.cell_to_ctype
   ctype_to_lvertex_to_qcoords = map(get_vertex_coordinates, polytopes)
@@ -290,10 +294,6 @@ function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation,glue)
   FaceCompressedVector(ctype_to_lface_to_pindex_to_qcoords,glue)
 end
 
-function _compute_face_to_q_vertex_coords(trian::BoundaryTriangulation{Dc,Dp,A,<:FaceToCellGlue}) where {Dc,Dp,A}
-  _compute_face_to_q_vertex_coords(trian,trian.glue)
-end
-
 struct FaceCompressedVector{T,G<:FaceToCellGlue} <: AbstractVector{T}
   ctype_lface_pindex_to_value::Vector{Vector{Vector{T}}}
   glue::G

src/Geometry/Triangulations.jl

@@ -52,10 +52,18 @@ get_cell_reffe(trian::Triangulation) = get_cell_reffe(get_grid(trian))
 is_first_order(trian::Triangulation) = is_first_order(get_grid(trian))
 
 # This is the most used glue, but others are possible, see e.g. SkeletonGlue.
-struct FaceToFaceGlue{A,B,C}
+# Dc is the topological dimension of the cells of the Triangulation{Dc} from
+# which the FaceToFaceGlue object is created
+struct FaceToFaceGlue{Dc,A,B,C}
   tface_to_mface::A
   tface_to_mface_map::B
   mface_to_tface::C
+  function FaceToFaceGlue(Dc::Integer,
+    tface_to_mface::A,
+    tface_to_mface_map::B,
+    mface_to_tface::C) where {A,B,C}
+    new{Dc,A,B,C}(tface_to_mface,tface_to_mface_map,mface_to_tface)
+  end
 end
 
 """
@@ -184,7 +192,7 @@ function get_glue(trian::BodyFittedTriangulation{Dt,Dp,A,B,C,Injective},::Val{Dt
     nmfaces = num_faces(trian.model,Dt)
     mface_to_tface = PosNegPartition(trian.tface_to_mface,Int32(nmfaces))
   end
-  FaceToFaceGlue(trian.tface_to_mface,tface_to_mface_map,mface_to_tface)
+  FaceToFaceGlue(Dt,trian.tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
 function get_glue(trian::BodyFittedTriangulation{Dt,Dp,A,B,C,NonInjective},::Val{Dt}) where {Dt,Dp,A,B,C}
@@ -193,7 +201,7 @@ function get_glue(trian::BodyFittedTriangulation{Dt,Dp,A,B,C,NonInjective},::Val
   # Whenever tface_to_mface is non-injective, we currently avoid the computation of 
   # mface_to_tface, which relies on PosNegPartition. This is a limitation that we should 
   # face in the future on those scenarios on which we need mface_to_tface.
-  FaceToFaceGlue(trian.tface_to_mface,tface_to_mface_map,mface_to_tface)
+  FaceToFaceGlue(Dt,trian.tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
 #function get_glue(trian::BodyFittedTriangulation{Dt},::Val{Dm}) where {Dt,Dm}
@@ -433,7 +441,7 @@ function _compose_glues(rglue,dglue)
   nothing
 end
 
-function _compose_glues(rglue::FaceToFaceGlue,dglue::FaceToFaceGlue)
+function _compose_glues(rglue::FaceToFaceGlue{Dc},dglue::FaceToFaceGlue) where Dc
   rface_to_mface = rglue.tface_to_mface
   dface_to_rface = dglue.tface_to_mface
   dface_to_mface = collect(lazy_map(Reindex(rface_to_mface),dface_to_rface))
@@ -442,7 +450,7 @@ function _compose_glues(rglue::FaceToFaceGlue,dglue::FaceToFaceGlue)
   dface_to_mface_map1 = lazy_map(Reindex(rface_to_mface_map),dface_to_rface)
   dface_to_mface_map = lazy_map(∘,dface_to_mface_map1,dface_to_rface_map)
   mface_to_dface = nothing
-  FaceToFaceGlue(dface_to_mface,dface_to_mface_map,mface_to_dface)
+  FaceToFaceGlue(Dc,dface_to_mface,dface_to_mface_map,mface_to_dface)
 end
 
 struct GenericTriangulation{Dc,Dp,A,B,C} <: Triangulation{Dc,Dp}
@@ -507,7 +515,7 @@ function get_glue(t::TriangulationView,::Val{d}) where d
   view(parent,t.cell_to_parent_cell)
 end
 
-function Base.view(glue::FaceToFaceGlue,ids::AbstractArray)
+function Base.view(glue::FaceToFaceGlue{Dc},ids::AbstractArray) where Dc
   tface_to_mface = lazy_map(Reindex(glue.tface_to_mface),ids)
   tface_to_mface_map = lazy_map(Reindex(glue.tface_to_mface_map),ids)
   if glue.mface_to_tface === nothing
@@ -516,7 +524,7 @@ function Base.view(glue::FaceToFaceGlue,ids::AbstractArray)
     nmfaces = length(glue.mface_to_tface)
     mface_to_tface = PosNegPartition(tface_to_mface,Int32(nmfaces))
   end
-  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
+  FaceToFaceGlue(Dc,tface_to_mface,tface_to_mface_map,mface_to_tface)
 end
 
 function get_facet_normal(trian::TriangulationView)

src/MultiField/MultiFieldFESpaces.jl

@@ -152,20 +152,21 @@ struct MultiFieldFEBasisComponent{B} <: FEBasis
   single_field::B
   fieldid::Int
   nfields::Int
-  function MultiFieldFEBasisComponent(
-    single_field::SingleFieldFEBasis,fieldid::Integer,nfields::Integer)
-    function block_dofs(cell_bs,::TestBasis,fieldid,nfields)
-      cell_basis = lazy_map(BlockMap(nfields,fieldid),cell_bs)
-    end
-    function block_dofs(cell_bs,::TrialBasis,fieldid,nfields)
-      cell_basis = lazy_map(BlockMap((1,nfields),fieldid),cell_bs)
-    end
-    B = typeof(single_field)
-    cell_bs = get_data(single_field)
-    bsty = BasisStyle(single_field)
-    cell_basis = block_dofs(cell_bs,bsty,fieldid,nfields)
-    new{B}(cell_basis,single_field,fieldid,nfields)
+end
+
+function MultiFieldFEBasisComponent(
+  single_field::SingleFieldFEBasis,fieldid::Integer,nfields::Integer)
+  function block_dofs(cell_bs,::TestBasis,fieldid,nfields)
+    cell_basis = lazy_map(BlockMap(nfields,fieldid),cell_bs)
+  end
+  function block_dofs(cell_bs,::TrialBasis,fieldid,nfields)
+    cell_basis = lazy_map(BlockMap((1,nfields),fieldid),cell_bs)
   end
+  B = typeof(single_field)
+  cell_bs = get_data(single_field)
+  bsty = BasisStyle(single_field)
+  cell_basis = block_dofs(cell_bs,bsty,fieldid,nfields)
+  MultiFieldFEBasisComponent{B}(cell_basis,single_field,fieldid,nfields)
 end
 
 CellData.get_data(f::MultiFieldFEBasisComponent) = f.cell_basis

src/ReferenceFEs/Quadratures.jl

@@ -99,13 +99,7 @@ get_weights(q::GenericQuadrature) = q.weights
 
 get_name(q::GenericQuadrature) = q.name
 
-function GenericQuadrature(a::Quadrature)
-  GenericQuadrature(get_coordinates(a),get_weights(a),get_name(a))
-end
 
-function GenericQuadrature(a::GenericQuadrature)
-  a
-end
 
 # Quadrature factory
 

test/ODEsTests/runtests.jl

@@ -18,7 +18,7 @@ using Test
 
 @time @testset "TransientFEOperatorsSolutions" begin include("TransientFEOperatorsSolutionsTests.jl") end
 
-@time @testset "TransientFEProblems" begin include("TransientFEProblemsTests.jl") end
+# @time @testset "TransientFEProblems" begin include("TransientFEProblemsTests.jl") end
 
 # @time @testset "DiffEqsWrappers" begin include("_DiffEqsWrappersTests.jl") end
 

test/runtests.jl

@@ -21,7 +21,7 @@ using Test
 @time @testset "Geometry" begin include("GeometryTests/runtests.jl") end
 
 @time @testset "CellData" begin include("CellDataTests/runtests.jl") end
-
+ 
 @time @testset "Visualization" begin include("VisualizationTests/runtests.jl") end
 
 @time @testset "FESpaces (1/2)" begin include("FESpacesTests/runtests_1.jl") end