Skip to content
Snippets Groups Projects
Select Git revision
  • e32f573f287ed41d8835dad9192519908054f0da
  • main default protected
  • askarpza-main-patch-76094
  • polynomial_regression
  • optimization
  • v0.8.0
  • v0.7.1
  • v0.7.0
  • v0.6.0
  • v0.5.0
  • v0.4.1
  • v0.4.0
  • v0.3.0
  • v0.2.0
14 results

test_grassmann.py

Blame
  • function.jl 2.07 KiB
    using Statistics: mean
    export Mapper, FeFunction, P1, DP0, DP1, interpolate!
    
    struct P1 end
    struct DP0 end
    struct DP1 end
    
    struct Mapper{T}
        mesh
        ndofs::Int
        map::Array{Int, 2}
    end
    
    struct FeFunction
        mapper::Mapper
        data::Vector{Float64}
        name::String
    end
    
    function FeFunction(mapper, name=string(gensym("f")))
        data = Vector{Float64}(undef, mapper.ndofs)
        return FeFunction(mapper, data, name)
    end
    
    # P1 Elements (1st order Lagrange)
    
    function Mapper{P1}(mesh)
        # special case for P1: dofs correspond to vertices
        ndofs = size(mesh.vertices, 2)
        map = copy(mesh.cells)
        return Mapper{P1}(mesh, ndofs, map)
    end
    
    # DP0 Elements (0th order discontinuous Lagrange)
    
    function Mapper{DP0}(mesh)
        # special case for DP0: dofs correspond to cells
        ndofs = size(mesh.cells, 2)
        map = collect(reshape(axes(mesh.cells, 2), 1, :))
        return Mapper{DP0}(mesh, ndofs, map)
    end
    
    interpolate!(dst::FeFunction, src::Function) = interpolate!(dst, dst.mapper, src)
    
    function interpolate!(dst::FeFunction, mapper::Mapper{P1}, src::Function)
        mesh = mapper.mesh
        for cid in axes(mesh.cells, 2)
    	for ldof in axes(mapper.map, 1)
    	    # TODO: how do ldofs generically map to position?
    	    xid = mesh.cells[ldof, cid]
    	    x = mesh.vertices[:, xid]
    
    	    # currently scalar
    	    fx = src(x)
    
    	    gdof = mapper.map[ldof, cid]
    	    dst.data[gdof] = fx
    	end
        end
    end
    
    function interpolate!(dst::FeFunction, mapper::Mapper{DP0}, src::Function)
        mesh = mapper.mesh
        for cid in axes(mesh.cells, 2)
    	vertices = mesh.vertices[:, mesh.cells[:, cid]]
    	centroid = mean(vertices, dims=2)
    
    	fx = src(centroid)
    
    	gdof = mapper.map[1, cid]
    	dst.data[gdof] = fx
        end
    end
    
    append_data!(vtkfile, f::FeFunction) = append_data!(vtkfile, f, f.mapper)
    
    function append_data!(vtkfile, f::FeFunction, ::Mapper{P1})
        # FIXME: vector-valued data
        fdata = reshape(f.data, 1, :)
        vtk_point_data(vtkfile, fdata, f.name)
    end
    
    function append_data!(vtkfile, f::FeFunction, ::Mapper{DP0})
        # FIXME: vector-valued data
        fdata = reshape(f.data, 1, :)
        vtk_cell_data(vtkfile, fdata, f.name)
    end