[NOT TO MERGE] Refactor Electro and ThermoElectro models

src/ComputationalModels/BoundaryConditions.jl

@@ -12,10 +12,6 @@ end
 
 getindex(bc::MultiFieldBC, i) = bc.BoundaryCondition[i]
 
-include("EvolutionFunctions.jl")
-include("CartesianTags.jl")
-include("FaceLabeling.jl")
-
 
 struct MultiFieldTC{A} <: TimedependentCondition
     vh::A # could be a multifield or single field

src/ComputationalModels/ComputationalModels.jl

@@ -24,6 +24,11 @@ using GridapGmsh: GmshDiscreteModel
 
 import Base.getindex
 
+# Deprecation exports
+using HyperFEM.DiscreteModeling
+Base.@deprecate_binding CartesianTags DiscreteModeling.CartesianTags
+Base.@deprecate_binding EvolutionFunctions DiscreteModeling.EvolutionFunctions
+
 include("BoundaryConditions.jl")
 export DirichletBC
 export NeumannBC

src/ComputationalModels/CartesianTags.jl → src/DiscreteModeling/CartesianTags.jl

@@ -2,6 +2,24 @@
 "Shortcuts for the tags of cartesian discrete models."
 module CartesianTags
 
+export face0YZ, face1YZ, faceX0Z, faceX1Z, faceXY0, faceXY1
+
+export edgeX00, edgeX10, edgeX01, edgeX11
+export edge0Y0, edge1Y0, edge0Y1, edge1Y1
+export edge00Z, edge10Z, edge01Z, edge11Z
+
+export corner000, corner100, corner010, corner110
+export corner001, corner101, corner011, corner111
+
+export edgeX00⁺, edgeX10⁺, edgeX01⁺, edgeX11⁺
+export edge0Y0⁺, edge1Y0⁺, edge0Y1⁺, edge1Y1⁺
+export edge00Z⁺, edge10Z⁺, edge01Z⁺, edge11Z⁺
+
+export face0YZ⁺, face1YZ⁺, faceX0Z⁺, faceX1Z⁺, faceXY0⁺, faceXY1⁺
+
+# Deprecated tags
+export faceX0, faceX1, faceY0, faceY1, faceZ0, faceZ1
+
 # --- Face tags ---
 
 "Tags indicating the face at plane X0."

src/DiscreteModeling/DiscreteModeling.jl

@@ -0,0 +1,20 @@
+
+"""
+A bundle of helper tools to work with discrete models in space/time.
+"""
+module DiscreteModeling
+
+using Gridap
+
+export CartesianTags
+export EvolutionFunctions
+export add_tag_from_vertex_filter!
+export aspect_ratio
+export element_size
+
+include("CartesianTags.jl")
+include("EvolutionFunctions.jl")
+include("FaceLabeling.jl")
+include("MeshDescriptor.jl")
+
+end

src/DiscreteModeling/FaceLabeling.jl

@@ -0,0 +1,28 @@
+
+"""
+Create a new tag from a geometry and a coordinate-based filter function.
+The filter function takes in vertex coordinates and returns a boolean values. A geometrical 
+entity is tagged if all its vertices pass the filter.
+
+# See also
+- `Gridap.Geometry.face_labeling_from_vertex_filter`
+- `Gridap.Geometry.merge!`
+"""
+function add_tag_from_vertex_filter!(
+            labels::Gridap.Geometry.FaceLabeling,
+            tag::String,
+            geometry::Gridap.Geometry.DiscreteModel,
+            filter::Function)
+  new_labels = Gridap.Geometry.face_labeling_from_vertex_filter(geometry.grid_topology, tag, filter)
+  merge!(labels, new_labels)
+end
+
+
+function add_tag_from_vertex_filter!(
+            labels::Gridap.Geometry.FaceLabeling,
+            geometry::Gridap.Geometry.DiscreteModel,
+            tag::String,
+            filter::Function)
+  @warn "This method is deprecated. Call the function with: (labels, tag, geometry, filter)"
+  add_tag_from_vertex_filter!(labels, tag, geometry, filter)
+end

src/DiscreteModeling/MeshDescriptor.jl

@@ -0,0 +1,88 @@
+
+"""
+Return the aspect ratio of the underlying cartesian elements as a string.
+This function is only available for an underlying `CartesianGrid`.
+
+# Example
+    aspect_ratio(Ω)             # "51:51:5"
+    aspect_ratio(Ω, tol=0.05)   # "10:10:1"
+    aspect_ratio(uh⁺, tol=0.1)  # "10:10:1" 
+"""
+function aspect_ratio(grid::CartesianGrid; tol=1e-6)
+  descriptor = Gridap.Geometry.get_cartesian_descriptor(grid)
+  sizes = descriptor.sizes
+
+  rel = sizes ./ minimum(sizes)
+  best = nothing
+  best_error = Inf
+  best_complexity = Inf
+
+  for d in 0:20
+    candidate = round.(Int, rel .* d)
+    any(candidate .== 0) && continue
+    approx = candidate ./ d
+    err = maximum(abs.(approx .- rel) ./ rel)
+    if err < best_error
+      complexity = max(candidate...)
+      if complexity < best_complexity
+        best = candidate
+        best_error = err
+        best_complexity = complexity
+      end
+    end
+  end
+
+  if best_error > tol  # exact fallback
+      denominators = denominator.(rationalize.(sizes))
+      least_mult = lcm(denominators...)
+      best = round.(Int, sizes .* least_mult)
+  end
+
+  join(best, ":")
+end
+
+function aspect_ratio(model::CartesianDiscreteModel; kwargs...)
+  aspect_ratio(get_grid(model); kwargs...)
+end
+
+function aspect_ratio(triangulation::Triangulation; kwargs...)
+  aspect_ratio(get_background_model(triangulation); kwargs...)
+end
+
+function aspect_ratio(f::CellField; kwargs...)
+  aspect_ratio(get_triangulation(f); kwargs...)
+end
+
+
+"""
+Return the element size for a cartesian mesh.
+This function is only available for an underlying `CartesianGrid`.
+
+# Example
+    element_size(model)   # Compute the diagonal
+    element_size(uh, :x)  # Get the x-size of the underlying grid
+"""
+function element_size(grid::CartesianGrid)
+  descriptor = Gridap.Geometry.get_cartesian_descriptor(grid)
+  sizes = descriptor.sizes
+  sqrt(sum(abs2, sizes))
+end
+
+function element_size(grid::CartesianGrid, direction)
+  descriptor = Gridap.Geometry.get_cartesian_descriptor(grid)
+  direction_indices = Dict(:x => 1, :y => 2, :z => 3)
+  index = direction_indices[direction]
+  descriptor.sizes[index]
+end
+
+function element_size(model::CartesianDiscreteModel, args...)
+  element_size(get_grid(model), args...)
+end
+
+function element_size(triangulation::Triangulation, args...)
+  element_size(get_background_model(triangulation), args...)
+end
+
+function element_size(f::CellField, args...)
+  element_size(get_triangulation(f), args...)
+end

src/Exports.jl

@@ -24,7 +24,6 @@ end
 @publish TensorAlgebra Tensorize
 
 
-@publish PhysicalModels DerivativeStrategy
 @publish PhysicalModels LinearElasticity3D
 @publish PhysicalModels LinearElasticity2D
 @publish PhysicalModels Yeoh3D
@@ -45,10 +44,12 @@ end
 @publish PhysicalModels NonlinearMooneyRivlin_CV
 @publish PhysicalModels NonlinearIncompressibleMooneyRivlin2D_CV
 @publish PhysicalModels EightChain
+@publish PhysicalModels EightChain5Terms
 @publish PhysicalModels TransverseIsotropy3D
 @publish PhysicalModels TransverseIsotropy2D
 @publish PhysicalModels ThermalModel
 @publish PhysicalModels ThermalVolumetric
+@publish PhysicalModels ThermalDeviatoric
 @publish PhysicalModels IdealDielectric
 @publish PhysicalModels Magnetic
 @publish PhysicalModels IdealMagnetic
@@ -58,6 +59,7 @@ end
 @publish PhysicalModels ElectroMechModel
 @publish PhysicalModels ThermoElectroMechModel
 @publish PhysicalModels ThermoMechModel
+@publish PhysicalModels ThermoElectroModel
 @publish PhysicalModels ThermoMech_Bonet
 @publish PhysicalModels ThermoMech_EntropicPolyconvex
 @publish PhysicalModels FlexoElectroModel
@@ -94,10 +96,11 @@ end
 @publish PhysicalModels getIsoInvariants
 
 @publish PhysicalModels ThermalLaw
+@publish PhysicalModels ConstantEnergyLaw
+@publish PhysicalModels ConstantCvLaw
 @publish PhysicalModels EntropicElasticityLaw
 @publish PhysicalModels NonlinearMeltingLaw
 @publish PhysicalModels NonlinearSofteningLaw
-@publish PhysicalModels TrigonometricLaw
 @publish PhysicalModels PolynomialLaw
 
 @publish PhysicalModels SecondPiola
@@ -152,8 +155,11 @@ end
 @publish ComputationalModels  TrialFESpace! # Exporting internal function of Gridap
 @publish ComputationalModels  L2_Projection
 
-# Note: the files FaceLabeling, CartesianTags and Evolution functions should be moved to a module different than ComputationalModels
-@publish ComputationalModels add_tag_from_vertex_filter!
+@publish DiscreteModeling CartesianTags
+@publish DiscreteModeling EvolutionFunctions
+@publish DiscreteModeling add_tag_from_vertex_filter!
+@publish DiscreteModeling aspect_ratio
+@publish DiscreteModeling element_size
 
 @publish Solvers IterativeSolver
 @publish Solvers Newton_RaphsonSolver

src/HyperFEM.jl

@@ -4,8 +4,8 @@ include("TensorAlgebra/TensorAlgebra.jl")
 include("PhysicalModels/PhysicalModels.jl")
 include("WeakForms/WeakForms.jl")
 include("Solvers/Solvers.jl")
+include("DiscreteModeling/DiscreteModeling.jl")
 include("ComputationalModels/ComputationalModels.jl")
-
 include("Io.jl")
 include("Exports.jl")
 

src/PhysicalModels/ElectricalModels.jl

@@ -9,3 +9,38 @@ struct IdealDielectric <: Electro
     new(ε)
   end
 end
+
+function (obj::IdealDielectric)()
+  J(F) = det(F)
+  H(F) = det(F) * inv(F)'
+
+  # Energy #
+  HE(F, E) = H(F) * E
+  HEHE(F, E) = HE(F, E) ⋅ HE(F, E)
+  Ψem(F, E) = (-obj.ε / (2.0 * J(F))) * HEHE(F, E)
+
+  # First Derivatives #
+  ∂Ψem_∂H(F, E) = (-obj.ε / (J(F))) * (HE(F, E) ⊗ E)
+  ∂Ψem_∂J(F, E) = (+obj.ε / (2.0 * J(F)^2.0)) * HEHE(F, E)
+  ∂Ψem_∂E(F, E) = (-obj.ε / (J(F))) * (H(F)' * HE(F, E))
+  ∂Ψem∂F(F, E) = ∂Ψem_∂H(F, E) × F + ∂Ψem_∂J(F, E) * H(F)
+  ∂Ψem∂E(F, E) = ∂Ψem_∂E(F, E)
+
+  # Second Derivatives #
+  ∂Ψem_HH(F, E) = (-obj.ε / (J(F))) * (I3 ⊗₁₃²⁴ (E ⊗ E))
+  ∂Ψem_HJ(F, E) = (+obj.ε / (J(F))^2.0) * (HE(F, E) ⊗ E)
+  ∂Ψem_JJ(F, E) = (-obj.ε / (J(F))^3.0) * HEHE(F, E)
+  ∂Ψem∂FF(F, E) = (F × (∂Ψem_HH(F, E) × F)) +
+                  H(F) ⊗₁₂³⁴ (∂Ψem_HJ(F, E) × F) +
+                  (∂Ψem_HJ(F, E) × F) ⊗₁₂³⁴ H(F) +
+                  ∂Ψem_JJ(F, E) * (H(F) ⊗₁₂³⁴ H(F)) +
+                  ×ᵢ⁴(∂Ψem_∂H(F, E) + ∂Ψem_∂J(F, E) * F)
+
+  ∂Ψem_EH(F, E) = (-obj.ε / (J(F))) * ((I3 ⊗₁₃² HE(F, E)) + (H(F)' ⊗₁₂³ E))
+  ∂Ψem_EJ(F, E) = (+obj.ε / (J(F))^2.0) * (H(F)' * HE(F, E))
+  ∂Ψem∂EF(F, E) = (∂Ψem_EH(F, E) × F) + (∂Ψem_EJ(F, E) ⊗₁²³ H(F))
+
+  ∂Ψem∂EE(F, E) = (-obj.ε / (J(F))) * (H(F)' * H(F))
+
+  return (Ψem, ∂Ψem∂F, ∂Ψem∂E, ∂Ψem∂FF, ∂Ψem∂EF, ∂Ψem∂EE)
+end