diff --git a/.gitignore b/.gitignore index 4fc850442d..3acf1cba2e 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,9 @@ *.jl.*.cov *.jl.mem *.vtu +*.vtk +*.pvtu +*.pvtk *.pvd examples/.ipynb_checkpoints/ docs/build/ @@ -10,3 +13,8 @@ docs/src/examples/ *.DS_Store /Manifest.toml /test/coverage/Manifest.toml +/benchmark/Manifest.toml + +LocalPreferences.toml +docs/LocalPreferences.toml +benchmark/tune.json diff --git a/Project.toml b/Project.toml index 928feae166..d41d176fb9 100644 --- a/Project.toml +++ b/Project.toml @@ -12,14 +12,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" -[weakdeps] -BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" -Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" - -[extensions] -FerriteBlockArrays = "BlockArrays" -FerriteMetis = "Metis" - [compat] BlockArrays = "0.16" EnumX = "1" @@ -31,6 +23,10 @@ Tensors = "1" WriteVTK = "1.13" julia = "1.6" +[extensions] +FerriteBlockArrays = "BlockArrays" +FerriteMetis = "Metis" + [extras] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" @@ -38,6 +34,7 @@ FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" @@ -46,7 +43,10 @@ SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" -Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" [targets] test = ["BlockArrays", "Downloads", "FerriteGmsh", "ForwardDiff", "Gmsh", "IterativeSolvers", "Metis", "NBInclude", "ProgressMeter", "Random", "SHA", "StableRNGs", "Test", "TimerOutputs", "Logging"] + +[weakdeps] +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" +Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" diff --git a/benchmark/Makefile b/benchmark/Makefile new file mode 100644 index 0000000000..96bf912a30 --- /dev/null +++ b/benchmark/Makefile @@ -0,0 +1,26 @@ +# Allow custom julia commands +JULIA_CMD ?= julia + +BENCHMARK_DIR := $(shell dirname $(abspath $(firstword $(MAKEFILE_LIST)))) + +default: benchmark + +clean: + rm Manifest.toml + +Manifest.toml: + ${JULIA_CMD} --project=${BENCHMARK_DIR} -e 'using Pkg; Pkg.develop(path=".."); Pkg.instantiate(); Pkg.precompile();' + +benchmark: Manifest.toml + ${JULIA_CMD} --project=${BENCHMARK_DIR} ${BENCHMARK_DIR}/runbenchmarks.jl benchmark/benchmarks.jl + +benchmark-%: Manifest.toml + ${JULIA_CMD} --project=${BENCHMARK_DIR} ${BENCHMARK_DIR}/runbenchmarks.jl benchmark/benchmarks.jl ${@:benchmark-%=%} + +compare: Manifest.toml + ${JULIA_CMD} --project=${BENCHMARK_DIR} ${BENCHMARK_DIR}/runcomparison.jl ${target} ${baseline} benchmark/benchmarks.jl + +compare-%: Manifest.toml + ${JULIA_CMD} --project=${BENCHMARK_DIR} ${BENCHMARK_DIR}/runcomparison.jl ${target} ${baseline} benchmark/benchmarks.jl ${@:benchmark-%=%} + +.PHONY: default clean benchmark benchmark-% compare compare-% diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000000..163643236e --- /dev/null +++ b/benchmark/Project.toml @@ -0,0 +1,4 @@ +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992" +PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d" diff --git a/benchmark/benchmarks-assembly.jl b/benchmark/benchmarks-assembly.jl new file mode 100644 index 0000000000..3c646a30db --- /dev/null +++ b/benchmark/benchmarks-assembly.jl @@ -0,0 +1,70 @@ +#----------------------------------------------------------------------# +# Assembly functionality benchmarks +#----------------------------------------------------------------------# +SUITE["assembly"] = BenchmarkGroup() + +# Permute over common combinations for some commonly required local matrices. +SUITE["assembly"]["common-local"] = BenchmarkGroup() +COMMON_LOCAL_ASSEMBLY = SUITE["assembly"]["common-local"] +for spatial_dim ∈ 1:3 + ξ_dummy = Vec{spatial_dim}(ntuple(x->0.0, spatial_dim)) + COMMON_LOCAL_ASSEMBLY["spatial-dim",spatial_dim] = BenchmarkGroup() + for geo_type ∈ FerriteBenchmarkHelper.geo_types_for_spatial_dim(spatial_dim) + COMMON_LOCAL_ASSEMBLY["spatial-dim",spatial_dim][string(geo_type)] = BenchmarkGroup() + + grid = generate_grid(geo_type, tuple(repeat([1], spatial_dim)...)); + ref_type = FerriteBenchmarkHelper.default_refshape(geo_type) + ip_geo = Ferrite.default_interpolation(geo_type) + + # Nodal interpolation tests + for order ∈ 1:2, ip_type ∈ [Lagrange, Serendipity] + ip = ip_type{spatial_dim, ref_type, order}() + + # Skip over elements which are not implemented + !applicable(Ferrite.value, ip, 1, ξ_dummy) && continue + + qr = QuadratureRule{spatial_dim, ref_type}(2*order-1) + + # Currenlty we just benchmark nodal Lagrange bases. + COMMON_LOCAL_ASSEMBLY["spatial-dim",spatial_dim][string(geo_type)][string(ip_type),string(order)] = BenchmarkGroup() + LAGRANGE_SUITE = COMMON_LOCAL_ASSEMBLY["spatial-dim",spatial_dim][string(geo_type)][string(ip_type),string(order)] + LAGRANGE_SUITE["fe-values"] = BenchmarkGroup() + LAGRANGE_SUITE["ritz-galerkin"] = BenchmarkGroup() + LAGRANGE_SUITE["petrov-galerkin"] = BenchmarkGroup() + + # Note: at the time of writing this PR the ctor makes the heavy lifting and caches important values. + LAGRANGE_SUITE["fe-values"]["scalar"] = @benchmarkable CellScalarValues($qr, $ip, $ip_geo); + LAGRANGE_SUITE["fe-values"]["vector"] = @benchmarkable CellVectorValues($qr, $ip, $ip_geo); + + csv = CellScalarValues(qr, ip, ip_geo); + csv2 = CellScalarValues(qr, ip, ip_geo); + + cvv = CellVectorValues(qr, ip, ip_geo); + cvv2 = CellVectorValues(qr, ip, ip_geo); + + # Scalar shape φ and test ψ: ∫ φ ψ + LAGRANGE_SUITE["ritz-galerkin"]["mass"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $csv, shape_value, shape_value, *) + LAGRANGE_SUITE["petrov-galerkin"]["mass"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $csv, shape_value, $csv2, shape_value, *) + # Vectorial shape φ and test ψ: ∫ φ ⋅ ψ + LAGRANGE_SUITE["ritz-galerkin"]["vector-mass"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $cvv, shape_value, shape_value, ⋅) + LAGRANGE_SUITE["petrov-galerkin"]["vector-mass"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $cvv, shape_value, $cvv2, shape_value, ⋅) + # Scalar shape φ and test ψ: ∫ ∇φ ⋅ ∇ψ + LAGRANGE_SUITE["ritz-galerkin"]["Laplace"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $csv, shape_gradient, shape_gradient, ⋅) + LAGRANGE_SUITE["petrov-galerkin"]["Laplace"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $csv, shape_gradient, $csv2, shape_gradient, ⋅) + # Vectorial shape φ and test ψ: ∫ ∇φ : ∇ψ + LAGRANGE_SUITE["ritz-galerkin"]["vector-Laplace"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $cvv, shape_gradient, shape_gradient, ⊡) + LAGRANGE_SUITE["petrov-galerkin"]["vector-Laplace"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $cvv, shape_gradient, $cvv2, shape_gradient, ⊡) + # Vectorial shape φ and scalar test ψ: ∫ (∇ ⋅ φ) ψ + LAGRANGE_SUITE["petrov-galerkin"]["pressure-velocity"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $cvv, shape_divergence, $csv, shape_value, *) + + if spatial_dim > 1 + qr_face = QuadratureRule{spatial_dim-1, ref_type}(2*order-1) + fsv = FaceScalarValues(qr_face, ip, ip_geo); + fsv2 = FaceScalarValues(qr_face, ip, ip_geo); + + LAGRANGE_SUITE["ritz-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, shape_value, *) + LAGRANGE_SUITE["petrov-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, $fsv2, shape_value, *) + end + end + end +end diff --git a/benchmark/benchmarks-boundary-conditions.jl b/benchmark/benchmarks-boundary-conditions.jl new file mode 100644 index 0000000000..b2a3b0ac7d --- /dev/null +++ b/benchmark/benchmarks-boundary-conditions.jl @@ -0,0 +1,44 @@ +#----------------------------------------------------------------------# +# Boundary condition benchmarks +#----------------------------------------------------------------------# +SUITE["boundary-conditions"] = BenchmarkGroup() + +SUITE["boundary-conditions"]["Dirichlet"] = BenchmarkGroup() +DIRICHLET_SUITE = SUITE["boundary-conditions"]["Dirichlet"] +# span artifical scope... +for spatial_dim ∈ [2] + # Benchmark application on global system + DIRICHLET_SUITE["global"] = BenchmarkGroup() + + geo_type = Quadrilateral + grid = generate_grid(geo_type, ntuple(x->2, spatial_dim)); + ref_type = FerriteBenchmarkHelper.default_refshape(geo_type) + ip_geo = Ferrite.default_interpolation(geo_type) + order = 2 + + # assemble a mass matrix to apply BCs on (because its cheap) + ip = Lagrange{spatial_dim, ref_type, order}() + qr = QuadratureRule{spatial_dim, ref_type}(2*order-1) + cellvalues = CellScalarValues(qr, ip, ip_geo); + dh = DofHandler(grid) + push!(dh, :u, 1, ip) + close!(dh); + + ch = ConstraintHandler(dh); + ∂Ω = union(getfaceset.((grid, ), ["left"])...); + dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0) + add!(ch, dbc); + close!(ch); + + # Non-symmetric application + M, f = FerriteAssemblyHelper._assemble_mass(dh, cellvalues, false); + DIRICHLET_SUITE["global"]["apply!(M,f,APPLY_TRANSPOSE)"] = @benchmarkable apply!($M, $f, $ch; strategy=$(Ferrite.APPLY_TRANSPOSE)); + DIRICHLET_SUITE["global"]["apply!(M,f,APPLY_INPLACE)"] = @benchmarkable apply!($M, $f, $ch; strategy=$(Ferrite.APPLY_INPLACE)); + # Symmetric application + M, f = FerriteAssemblyHelper._assemble_mass(dh, cellvalues, true); + DIRICHLET_SUITE["global"]["apply!(M_sym,f,APPLY_TRANSPOSE)"] = @benchmarkable apply!($M, $f, $ch; strategy=$(Ferrite.APPLY_TRANSPOSE)); + DIRICHLET_SUITE["global"]["apply!(M_sym,f,APPLY_INPLACE)"] = @benchmarkable apply!($M, $f, $ch; strategy=$(Ferrite.APPLY_INPLACE)); + + DIRICHLET_SUITE["global"]["apply!(f)"] = @benchmarkable apply!($f, $ch); + DIRICHLET_SUITE["global"]["apply_zero!(f)"] = @benchmarkable apply!($f, $ch); +end diff --git a/benchmark/benchmarks-dofs.jl b/benchmark/benchmarks-dofs.jl new file mode 100644 index 0000000000..33106120ee --- /dev/null +++ b/benchmark/benchmarks-dofs.jl @@ -0,0 +1,92 @@ + +#----------------------------------------------------------------------# +# Benchmarks around the dof management +#----------------------------------------------------------------------# +SUITE["dof-management"] = BenchmarkGroup() +SUITE["dof-management"]["numbering"] = BenchmarkGroup() +# !!! NOTE close! must wrapped into a custom function, because consecutive calls to close!, since the dofs are already distributed. +NUMBERING_SUITE = SUITE["dof-management"]["numbering"] +for spatial_dim ∈ [3]# 1:3 + NUMBERING_SUITE["spatial-dim",spatial_dim] = BenchmarkGroup() + for geo_type ∈ FerriteBenchmarkHelper.geo_types_for_spatial_dim(spatial_dim) + NUMBERING_SUITE["spatial-dim",spatial_dim][string(geo_type)] = BenchmarkGroup() + + ref_type = FerriteBenchmarkHelper.default_refshape(geo_type) + + for grid_size ∈ [2]#[3, 6, 9] #multiple grid sized to estimate computational complexity... + NUMBERING_SUITE["spatial-dim",spatial_dim][string(geo_type)]["grid-size-",grid_size] = BenchmarkGroup() + NUMBERING_SUITE["spatial-dim",spatial_dim][string(geo_type)]["grid-size-",grid_size] = BenchmarkGroup() + + grid = generate_grid(geo_type, ntuple(x->grid_size, spatial_dim)); + + for field_dim ∈ [3]#1:3 + NUMBERING_SUITE["spatial-dim",spatial_dim][string(geo_type)]["grid-size-",grid_size]["field-dim-", field_dim] = BenchmarkGroup() + NUMBERING_FIELD_DIM_SUITE = NUMBERING_SUITE["spatial-dim",spatial_dim][string(geo_type)]["grid-size-",grid_size]["field-dim-", field_dim] + # Lagrange tests + for order ∈ 1:2 + ip = Lagrange{spatial_dim, ref_type, order}() + + # Skip over elements which are not implemented + ξ_dummy = Vec{spatial_dim}(ntuple(x->0.0, spatial_dim)) + !applicable(Ferrite.value, ip, 1, ξ_dummy) && continue + !applicable(Ferrite.value, ip, 1, ξ_dummy) && continue + + NUMBERING_FIELD_DIM_SUITE["Lagrange",order] = BenchmarkGroup() + LAGRANGE_SUITE = NUMBERING_FIELD_DIM_SUITE["Lagrange",order] + order2 = max(order-1, 1) + ip2 = Lagrange{spatial_dim, ref_type, order2}() + + LAGRANGE_SUITE["DofHandler"] = BenchmarkGroup() + + close_helper = function(grid, ip) + dh = DofHandler(grid) + push!(dh, :u, field_dim, ip) + close!(dh) + end + LAGRANGE_SUITE["DofHandler"]["one-field"] = @benchmarkable $close_helper($grid, $ip) + + close_helper = function(grid, ip, ip2) + dh = DofHandler(grid) + push!(dh, :u, field_dim, ip) + push!(dh, :p, 1, ip2) + close!(dh) + end + LAGRANGE_SUITE["DofHandler"]["two-fields"] = @benchmarkable $close_helper($grid, $ip, $ip2) + + + LAGRANGE_SUITE["MixedDofHandler"] = BenchmarkGroup() + f1 = Field(:u, ip, field_dim) + f2 = Field(:p, ip2, 1) + + close_helper = function(grid, f1) + dh = MixedDofHandler(grid) + push!(dh, FieldHandler([f1], Set(1:getncells(grid)))) + close!(dh) + end + LAGRANGE_SUITE["MixedDofHandler"]["one-field"] = @benchmarkable $close_helper($grid, $f1) + + close_helper = function(grid, f1) + dh = MixedDofHandler(grid) + push!(dh, FieldHandler([f1], Set(1:Int(round(getncells(grid)/2))))) + close!(dh) + end + LAGRANGE_SUITE["MixedDofHandler"]["one-field-subdomain"] = @benchmarkable $close_helper($grid, $f1) + + close_helper = function(grid, f1, f2) + dh = MixedDofHandler(grid) + push!(dh, FieldHandler([f1, f2], Set(1:getncells(grid)))) + close!(dh) + end + LAGRANGE_SUITE["MixedDofHandler"]["two-fields"] = @benchmarkable $close_helper($grid, $f1, $f2) + + close_helper = function(grid, f1, f2) + dh = MixedDofHandler(grid) + push!(dh, FieldHandler([f1, f2], Set(1:Int(round(getncells(grid)/2))))) + close!(dh) + end + LAGRANGE_SUITE["MixedDofHandler"]["two-fields-subdomain"] = @benchmarkable $close_helper($grid, $f1, $f2) + end + end + end + end +end diff --git a/benchmark/benchmarks-mesh.jl b/benchmark/benchmarks-mesh.jl new file mode 100644 index 0000000000..34d7413f38 --- /dev/null +++ b/benchmark/benchmarks-mesh.jl @@ -0,0 +1,21 @@ + +#----------------------------------------------------------------------# +# Benchmarks for mesh functionality within Ferrite +#----------------------------------------------------------------------# +SUITE["mesh"] = BenchmarkGroup() + +# Generator benchmarks +SUITE["mesh"]["generator"] = BenchmarkGroup() + +# Strucutred hyperrectangle generators +SUITE["mesh"]["generator"]["hyperrectangle"] = BenchmarkGroup() +HYPERRECTANGLE_GENERATOR = SUITE["mesh"]["generator"]["hyperrectangle"] +for spatial_dim ∈ 1:3 + HYPERRECTANGLE_GENERATOR["spatial-dim",spatial_dim] = BenchmarkGroup() + for geo_type ∈ FerriteBenchmarkHelper.geo_types_for_spatial_dim(spatial_dim) + HYPERRECTANGLE_GENERATOR["spatial-dim",spatial_dim][string(geo_type)] = @benchmarkable generate_grid($geo_type, $(ntuple(x->4, spatial_dim))); + end +end + +# TODO AMR performance +# TODO topology performance diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 0000000000..683da055f8 --- /dev/null +++ b/benchmark/benchmarks.jl @@ -0,0 +1,22 @@ +using BenchmarkTools +using Ferrite + +const selected = get(ENV, "FERRITE_SELECTED_BENCHMARKS", "all") +const runall = selected == "all" + +include("helper.jl") + +const SUITE = BenchmarkGroup() + +if runall || selected == "mesh" + include("benchmarks-mesh.jl") +end +if runall || selected == "dofs" + include("benchmarks-dofs.jl") +end +if runall || selected == "assembly" + include("benchmarks-assembly.jl") +end +if runall || selected == "boundary-conditions" + include("benchmarks-boundary-conditions.jl") +end diff --git a/benchmark/helper.jl b/benchmark/helper.jl new file mode 100644 index 0000000000..b3f8d7d1e6 --- /dev/null +++ b/benchmark/helper.jl @@ -0,0 +1,159 @@ +module FerriteBenchmarkHelper + +using Ferrite + +function geo_types_for_spatial_dim(spatial_dim) + spatial_dim == 1 && return [Line, QuadraticLine] + spatial_dim == 2 && return [Triangle, QuadraticTriangle, Quadrilateral, QuadraticQuadrilateral] + spatial_dim == 3 && return [Tetrahedron, Hexahedron] # Quadratic* not yet functional in 3D. 3D triangle missing. Embedded also missing. +end + +default_refshape(t::Type{C}) where {C <: Ferrite.AbstractCell} = typeof(Ferrite.default_interpolation(t)).parameters[2] + +end + + +module FerriteAssemblyHelper + +using Ferrite + +# Minimal Ritz-Galerkin type local assembly loop. +function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues::CellValues{dim,T,refshape}, f_shape, f_test, op) where {dim,T,refshape} + n_basefuncs = getnbasefunctions(cellvalues) + + Ke = zeros(n_basefuncs, n_basefuncs) + + X = getcoordinates(grid, 1) + reinit!(cellvalues, X) + + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + for i in 1:n_basefuncs + test = f_test(cellvalues, q_point, i) + for j in 1:n_basefuncs + shape = f_shape(cellvalues, q_point, j) + Ke[i, j] += op(test, shape) * dΩ + end + end + end + + Ke +end + +function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, facevalues::FaceValues{dim,T,refshape}, f_shape, f_test, op) where {dim,T,refshape} + n_basefuncs = getnbasefunctions(facevalues) + + f = zeros(n_basefuncs) + + X = getcoordinates(grid, 1) + for face in 1:nfaces(getcells(grid)[1]) + reinit!(facevalues, X, face) + + for q_point in 1:getnquadpoints(facevalues) + n = getnormal(facevalues, q_point) + dΓ = getdetJdV(facevalues, q_point) + for i in 1:n_basefuncs + test = f_test(facevalues, q_point, i) + for j in 1:n_basefuncs + shape = f_shape(facevalues, q_point, j) + f[i] += op(test, shape) ⋅ n * dΓ + end + end + end + end + + f +end + +# Minimal Petrov-Galerkin type local assembly loop. We assume that both function spaces share the same integration rule. Test is applied from the left. +function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues_shape::CellValues{dim,T,refshape}, f_shape, cellvalues_test::CellValues{dim,T,refshape}, f_test, op) where {dim,T,refshape} + n_basefuncs_shape = getnbasefunctions(cellvalues_shape) + n_basefuncs_test = getnbasefunctions(cellvalues_test) + + Ke = zeros(n_basefuncs_test, n_basefuncs_shape) + + #implicit assumption: Same geometry! + X_shape = zeros(Vec{dim,Float64}, Ferrite.getngeobasefunctions(cellvalues_shape)) + getcoordinates!(X_shape, grid, 1) + reinit!(cellvalues_shape, X_shape) + + X_test = zeros(Vec{dim,Float64}, Ferrite.getngeobasefunctions(cellvalues_test)) + getcoordinates!(X_test, grid, 1) + reinit!(cellvalues_test, X_test) + + for q_point in 1:getnquadpoints(cellvalues_test) #assume same quadrature rule + dΩ = getdetJdV(cellvalues_test, q_point) + for i in 1:n_basefuncs_test + test = f_test(cellvalues_test, q_point, i) + for j in 1:n_basefuncs_shape + shape = f_shape(cellvalues_shape, q_point, j) + Ke[i, j] += op(test, shape) * dΩ + end + end + end + + Ke +end + +function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, facevalues_shape::FaceValues{dim,T,refshape}, f_shape, facevalues_test::FaceValues{dim,T,refshape}, f_test, op) where {dim,T,refshape} + n_basefuncs_shape = getnbasefunctions(facevalues_shape) + n_basefuncs_test = getnbasefunctions(facevalues_test) + + f = zeros(n_basefuncs_test) + + X_shape = getcoordinates(grid, 1) + X_test = getcoordinates(grid, 1) + for face in 1:nfaces(getcells(grid)[1]) + reinit!(facevalues_shape, X_shape, face) + reinit!(facevalues_test, X_test, face) + + for q_point in 1:getnquadpoints(facevalues_shape) + n = getnormal(facevalues_test, q_point) + dΓ = getdetJdV(facevalues_test, q_point) + for i in 1:n_basefuncs_test + test = f_test(facevalues_test, q_point, i) + for j in 1:n_basefuncs_shape + shape = f_shape(facevalues_shape, q_point, j) + f[i] += op(test, shape) ⋅ n * dΓ + end + end + end + end + + f +end + +function _assemble_mass(dh, cellvalues, sym) + n_basefuncs = getnbasefunctions(cellvalues) + Me = zeros(n_basefuncs, n_basefuncs) + fe = zeros(n_basefuncs) + + M = sym ? create_symmetric_sparsity_pattern(dh) : create_sparsity_pattern(dh); + f = zeros(ndofs(dh)) + + assembler = start_assemble(M, f); + @inbounds for cell in CellIterator(dh) + fill!(Me, 0) + + reinit!(cellvalues, cell) + + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + + for i in 1:n_basefuncs + φ = shape_value(cellvalues, q_point, i) + fe[i] += φ * dΩ + for j in 1:n_basefuncs + ψ = shape_value(cellvalues, q_point, j) + Me[i, j] += (φ * ψ) * dΩ + end + end + end + + assemble!(assembler, celldofs(cell), Me, fe) + end + + return M, f +end + +end diff --git a/benchmark/runbenchmarks.jl b/benchmark/runbenchmarks.jl new file mode 100644 index 0000000000..d0e19cf627 --- /dev/null +++ b/benchmark/runbenchmarks.jl @@ -0,0 +1,21 @@ +if !(length(ARGS) == 1 || length(ARGS) == 2) + error("Usage: runbenchmarks.jl