Source code for openpnm.io._comsol
import numpy as np
import openpnm as op
[docs]
def network_to_comsol(network, filename=None):
r"""
Saves the network and geometry data from the given objects into the
specified file. This exports in 2D only where throats and pores
have rectangular and circular shapes, respectively.
Parameters
----------
network : Network
The network containing the desired data
Notes
-----
- The exported files contain COMSOL geometry objects, not meshes.
- This class exports in 2D only.
"""
if op.topotools.dimensionality(network).sum() == 3:
raise Exception("COMSOL I/O class only works for 2D networks!")
filename = network.name if filename is None else filename
f = open(f"{filename}.mphtxt", 'w')
header(file=f, Nr=network.Nt, Nc=network.Np)
cn = network['throat.conns']
p1 = network['pore.coords'][cn[:, 0]]
p2 = network['pore.coords'][cn[:, 1]]
# Compute the rotation angle of throats
dif_x = p2[:, 0]-p1[:, 0]
dif_y = p2[:, 1]-p1[:, 1]
# Avoid division by 0
m = np.array([dif_x_i != 0 for dif_x_i in dif_x])
r = np.zeros((len(dif_x)))
r[~m] = np.inf
r[m] = dif_y[m]/dif_x[m]
angles = np.arctan(r)
r_w = network['throat.diameter']
rectangles(file=f, pores1=p1, pores2=p2, alphas=angles, widths=r_w)
c_c = network['pore.coords']
c_r = network['pore.diameter']/2.0
circles(file=f, centers=c_c, radii=c_r)
f.close()
def header(file, Nr, Nc):
f = file
f.write('# Geometry exported by OpenPNM'+2*'\n')
f.write('# Major & minor version'+'\n')
f.write('0 1'+2*'\n')
f.write(str(Nc+Nr)+' '+'# number of tags'+'\n')
f.write('# Tags'+'\n')
for r in range(1, Nr+1):
tag = 'r'+str(int(r))
tag = str(len(tag))+' '+tag
f.write(tag+'\n')
for c in range(1, Nc+1):
tag = 'c'+str(int(c))
tag = str(len(tag))+' '+tag
f.write(tag+'\n')
f.write('\n'+str(Nc+Nr)+' '+'# number of types'+'\n')
f.write('# Types'+'\n')
for i in range(Nc+Nr):
f.write('3 obj'+'\n')
f.write('\n')
def rectangles(file, pores1, pores2, alphas, widths):
f = file
p1x = pores1[:, 0] + (widths/2)*np.sin(alphas)
p1y = pores1[:, 1] - (widths/2)*np.cos(alphas)
p2x = pores2[:, 0] + (widths/2)*np.sin(alphas)
p2y = pores2[:, 1] - (widths/2)*np.cos(alphas)
p3x = pores2[:, 0] - (widths/2)*np.sin(alphas)
p3y = pores2[:, 1] + (widths/2)*np.cos(alphas)
p4x = pores1[:, 0] - (widths/2)*np.sin(alphas)
p4y = pores1[:, 1] + (widths/2)*np.cos(alphas)
for r in range(len(pores1)):
f.write('# --------- rectangle nbr '+str(r+1)+' ---------'+2*'\n')
f.write('0 0 1'+'\n')
f.write('5 Geom2 # class'+'\n')
f.write('2 # version'+'\n')
f.write('2 # type'+'\n')
f.write('1 # voidsLabeled'+'\n')
f.write('1e-010 # gtol'+'\n')
f.write('0.0001 # resTol'+2*'\n')
f.write('4 # number of vertices'+'\n')
f.write('# Vertices'+'\n')
f.write('# X Y dom tol'+'\n')
f.write(str(p1x[r])+' '+str(p1y[r])+' '+'-1 NAN'+'\n')
f.write(str(p2x[r])+' '+str(p2y[r])+' '+'-1 NAN'+'\n')
f.write(str(p3x[r])+' '+str(p3y[r])+' '+'-1 NAN'+'\n')
f.write(str(p4x[r])+' '+str(p4y[r])+' '+'-1 NAN'+2*'\n')
f.write('4 # number of edges'+'\n')
f.write('# Edges'+'\n')
f.write('# vtx1 vtx2 s1 s2 up down mfd tol'+'\n')
f.write('2 1 0 1 0 1 1 NAN'+'\n')
f.write('3 2 0 1 0 1 2 NAN'+'\n')
f.write('4 3 0 1 0 1 3 NAN'+'\n')
f.write('1 4 0 1 0 1 4 NAN'+2*'\n')
f.write('4 # number of manifolds'+'\n')
f.write('# Manifolds'+2*'\n')
f.write('# Manifold #0'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('1 0 # degrees'+'\n')
f.write('2 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
f.write(str(p2x[r])+' '+str(p2y[r])+' '+'1'+'\n')
f.write(str(p1x[r])+' '+str(p1y[r])+' '+'1'+2*'\n')
f.write('# Manifold #1'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('1 0 # degrees'+'\n')
f.write('2 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
f.write(str(p3x[r])+' '+str(p3y[r])+' '+'1'+'\n')
f.write(str(p2x[r])+' '+str(p2y[r])+' '+'1'+2*'\n')
f.write('# Manifold #2'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('1 0 # degrees'+'\n')
f.write('2 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
f.write(str(p4x[r])+' '+str(p4y[r])+' '+'1'+'\n')
f.write(str(p3x[r])+' '+str(p3y[r])+' '+'1'+2*'\n')
f.write('# Manifold #3'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('1 0 # degrees'+'\n')
f.write('2 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
f.write(str(p1x[r])+' '+str(p1y[r])+' '+'1'+'\n')
f.write(str(p4x[r])+' '+str(p4y[r])+' '+'1'+2*'\n')
f.write('# Attributes'+'\n')
f.write('0 # nof attributes'+2*'\n')
def circles(file, centers, radii):
f = file
p1x = centers[:, 0]-radii
p1y = centers[:, 1]
p2x = centers[:, 0]
p2y = centers[:, 1]-radii
p3x = centers[:, 0]+radii
p3y = centers[:, 1]
p4x = centers[:, 0]
p4y = centers[:, 1]+radii
for c in range(len(centers)):
f.write('# --------- circle nbr '+str(c+1)+' ---------'+2*'\n')
f.write('0 0 1'+'\n')
f.write('5 Geom2 # class'+'\n')
f.write('2 # version'+'\n')
f.write('2 # type'+'\n')
f.write('1 # voidsLabeled'+'\n')
f.write('1e-010 # gtol'+'\n')
f.write('0.0001 # resTol'+2*'\n')
f.write('4 # number of vertices'+'\n')
f.write('# Vertices'+'\n')
f.write('# X Y dom tol'+'\n')
# extreme left point of the circle
f.write(str(p1x[c])+' '+str(p1y[c])+' '+'-1 NAN'+'\n')
# extreme bottom point of the circle
f.write(str(p2x[c])+' '+str(p2y[c])+' '+'-1 NAN'+'\n')
# extreme right point of the circle
f.write(str(p3x[c])+' '+str(p3y[c])+' '+'-1 NAN'+'\n')
# extreme top point of the circle
f.write(str(p4x[c])+' '+str(p4y[c])+' '+'-1 NAN'+2*'\n')
f.write('4 # number of edges'+'\n')
f.write('# Edges'+'\n')
f.write('# vtx1 vtx2 s1 s2 up down mfd tol'+'\n')
f.write('2 1 0 1 0 1 1 NAN'+'\n')
f.write('3 2 0 1 0 1 2 NAN'+'\n')
f.write('4 3 0 1 0 1 3 NAN'+'\n')
f.write('1 4 0 1 0 1 4 NAN'+2*'\n')
f.write('4 # number of manifolds'+'\n')
f.write('# Manifolds'+2*'\n')
# bottom left quart of the circle
f.write('# Manifold #0'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('2 0 # degrees'+'\n')
f.write('3 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
# extreme bottom point of the circle
f.write(str(p2x[c])+' '+str(p2y[c])+' '+'1'+'\n')
f.write(str(p1x[c])+' '+str(p2y[c])+' '+'0.70710678118654746'+'\n')
# extreme left point of the circle
f.write(str(p1x[c])+' '+str(p1y[c])+' '+'1'+2*'\n')
# bottom right quart of the circle
f.write('# Manifold #1'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('2 0 # degrees'+'\n')
f.write('3 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
# extreme right point of the circle
f.write(str(p3x[c])+' '+str(p3y[c])+' '+'1'+'\n')
f.write(str(p3x[c])+' '+str(p2y[c])+' '+'0.70710678118654746'+'\n')
# extreme bottom point of the circle
f.write(str(p2x[c])+' '+str(p2y[c])+' '+'1'+2*'\n')
# top right quart of the circle
f.write('# Manifold #2'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('2 0 # degrees'+'\n')
f.write('3 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
# extreme top point of the circle
f.write(str(p4x[c])+' '+str(p4y[c])+' '+'1'+'\n')
f.write(str(p3x[c])+' '+str(p4y[c])+' '+'0.70710678118654746'+'\n')
# extreme right point of the circle
f.write(str(p3x[c])+' '+str(p3y[c])+' '+'1'+2*'\n')
# top left quart of the circle
f.write('# Manifold #3'+'\n')
f.write('11 BezierCurve # class'+'\n')
f.write('0 0 # version'+'\n')
f.write('2 # sdim'+'\n')
f.write('0 2 1 # transformation'+'\n')
f.write('2 0 # degrees'+'\n')
f.write('3 # number of control points'+'\n')
f.write('# control point coords and weights'+'\n')
# extreme left point of the circle
f.write(str(p1x[c])+' '+str(p1y[c])+' '+'1'+'\n')
f.write(str(p1x[c])+' '+str(p4y[c])+' '+'0.70710678118654746'+'\n')
# extreme top point of the circle
f.write(str(p4x[c])+' '+str(p4y[c])+' '+'1'+2*'\n')
f.write('# Attributes'+'\n')
f.write('0 # nof attributes'+2*'\n')
