primitive version of xdmf writing

src/JuliaFEM.jl

@@ -6,6 +6,7 @@ VERSION < v"0.4-" && using Docile
 
 # import solvers
 include("elasticity_solver.jl")
+include("xdmf.jl")
 #using elasticity_solver
 
 export Model, new_model, new_field, get_field, add_nodes, get_nodes, add_elements, get_elements

src/xdmf.jl

@@ -0,0 +1,105 @@
+# This file is a part of JuliaFEM. License is MIT: https://github.com/ovainola/JuliaFEM/blob/master/README.md
+module xdmf
+
+using Logging
+@Logging.configure(level=INFO)
+
+using LightXML
+
+VERSION < v"0.4-" && using Docile
+
+# i add docstrings later
+
+function xdmf_new_model(xdmf_version="2.1")
+    xdoc = XMLDocument()
+    xroot = create_root(xdoc, "Xdmf")
+    set_attribute(xroot, "xmlns:xi", "http://www.w3.org/2001/XInclude")
+    set_attribute(xroot, "Version", xdmf_version)
+    domain = new_child(xroot, "Domain")
+    return xdoc, domain
+end
+
+function xdmf_new_temporal_collection(model)
+    temporal_collection = new_child(model, "Grid")
+    set_attribute(temporal_collection, "CollectionType", "Temporal")
+    set_attribute(temporal_collection, "GridType", "Collection")
+    set_attribute(temporal_collection, "Name", "Collection")
+    geometry = new_child(temporal_collection, "Geometry")
+    set_attribute(geometry, "Type", "None")
+    topology = new_child(temporal_collection, "Topology")
+    set_attribute(topology, "Dimensions", "0")
+    set_attribute(topology, "Type", "NoTopology")
+    return temporal_collection
+end
+
+function xdmf_new_grid(temporal_collection; time=0)
+    grid = new_child(temporal_collection, "Grid")
+    set_attribute(grid, "Name", "Grid")
+    time_ = new_child(grid, "Time")
+    set_attribute(time_, "Value", time)
+    return grid
+end
+
+function xdmf_new_mesh(grid, X, elmap)
+    geometry = new_child(grid, "Geometry")
+    set_attribute(geometry, "Type", "XYZ")
+    dataitem = new_child(geometry, "DataItem")
+    set_attribute(dataitem, "DataType", "Float")
+    set_attribute(dataitem, "Dimensions", length(X))
+    set_attribute(dataitem, "Format", "XML")
+    set_attribute(dataitem, "Precision", "4")
+    add_text(dataitem, join(X, " "))
+
+    topology = new_child(grid, "Topology")
+    set_attribute(topology, "Dimensions", "1")
+    set_attribute(topology, "Type", "Mixed")
+    dataitem = new_child(topology, "DataItem")
+    set_attribute(dataitem, "DataType", "Int")
+    set_attribute(dataitem, "Dimensions", length(elmap))
+    set_attribute(dataitem, "Format", "XML")
+    set_attribute(dataitem, "Precision", 4)
+    elmap2 = copy(elmap)
+    elmap2[2:end,:] -= 1
+    add_text(dataitem, join(elmap2, " "))
+end
+
+function xdmf_new_field(grid, name, source, data)
+    loc = Dict("elements" => "Cell",
+               "nodes" => "Node")
+
+    dim1, dim2 = size(data'')
+    if dim1 == 1
+        Type = "Scalar"
+    end
+    if dim1 == 3
+        Type = "Vector"
+    end
+
+    typ = string(typeof(data))
+    datatype = "unknown"
+    @debug("typeof: ", typ)
+    for j in ["Int", "Float"]
+        @debug(j)
+        if contains(typ, j)
+            datatype = j
+        end
+    end
+    if datatype == "unknown"
+        throw("unknown data type ", typ)
+    end
+
+
+    attribute = new_child(grid, "Attribute")
+    set_attribute(attribute, "Center", loc[source])
+    set_attribute(attribute, "Name", name)
+    set_attribute(attribute, "Type", Type)
+    dataitem = new_child(attribute, "DataItem")
+    set_attribute(dataitem, "DataType", datatype)
+    set_attribute(dataitem, "Dimensions", length(data))
+    set_attribute(dataitem, "Format", "XML")
+    set_attribute(dataitem, "Precision", 4)
+    add_text(dataitem, join(data, " "))
+end
+
+
+end

test/test_xdmf.jl

@@ -23,7 +23,7 @@ testdata = """\
           <DataItem DataType="Int" Dimensions="1" Format="XML" Precision="4">56</DataItem>
         </Attribute>
         <Attribute Center="Node" Name="Density" Type="Vector">
-          <DataItem DataType="Float" Dimensions="4" Format="XML" Precision="4">1 2 4 7</DataItem>
+          <DataItem DataType="Float" Dimensions="12" Format="XML" Precision="4">1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0</DataItem>
         </Attribute>
       </Grid>
     </Grid>
@@ -31,6 +31,9 @@ testdata = """\
 </Xdmf>
 """
 
+# makes no sense for vector field
+# <DataItem DataType="Float" Dimensions="4" Format="XML" Precision="4">1 2 4 7</DataItem>
+
 facts("test test data") do
   xdoc = XMLDocument()
   xroot = create_root(xdoc, "Xdmf")
@@ -85,10 +88,26 @@ facts("test test data") do
   set_attribute(attribute, "Type", "Vector")
   dataitem = new_child(attribute, "DataItem")
   set_attribute(dataitem, "DataType", "Float")
-  set_attribute(dataitem, "Dimensions", 4)
+  set_attribute(dataitem, "Dimensions", 12)
   set_attribute(dataitem, "Format", "XML")
   set_attribute(dataitem, "Precision", 4)
-  add_text(dataitem, "1 2 4 7")
+  add_text(dataitem, "1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0")
   #save_file(xdoc, "/tmp/model.xmf")
   @fact string(xdoc) => testdata
 end
+
+using JuliaFEM.xdmf: xdmf_new_model, xdmf_new_grid, xdmf_new_temporal_collection, xdmf_new_mesh, xdmf_new_field
+
+facts("test model write XML") do
+  X = [0 0 0; 1 0 0; 0 1 0; 0 0 1]'
+  elmap = [6 1 2 3 4]' # element type 6, nodes 1 2 3 4
+  temperature_field = [56]
+  density_field = Float64[1 2 3; 4 5 6; 7 8 9; 10 11 12]'
+  xdoc, model = xdmf_new_model()
+  temporal_collection = xdmf_new_temporal_collection(model)
+  grid = xdmf_new_grid(temporal_collection; time=123)
+  xdmf_new_mesh(grid, X, elmap)
+  xdmf_new_field(grid, "Temperature field", "elements", temperature_field)
+  xdmf_new_field(grid, "Density", "nodes", density_field)
+  @fact string(xdoc) => testdata
+end