Material anisotropy along a curve?

[JKM]

Hello,
For an electro ohmic quasi-static simulation, there is a check mark to specify a material setting as anisotropic, which will enable fields to enter conductivity on and off the diagonal of xx/yy/zz. As far as I understand it, this works well for a geometry with a fixed orientation/axis, but is there a way to specify the anisotropy with respect to a spline used to create the geometry (such as for muscles or nerves)? One might imagine this as using a local coordinate system rather than the global coordinate system.
Thank you for the help

[tcs]

I have a similar requirement at the moment. Did you get anywhere with it?

[bryn]

For this you would need to use the homogeneous anisotropic conductivity. The easiest way to use that feature I think, may be to create the field distribution in Python as a FloatFieldData on a RectilinearGrid and add this to the Analysis tree via the "TrivialProducer". You can then select the anisotropic conductivity from the Analysis pipeline.

Anisotropic conductivity distributions along a spline are sometimes computed by solving a Poisson equation in a tube, with Dirichlet boundary conditions at the top and bottom. This gives a main direction for the anisotropic tensor along the spline. You can create a tube from a spline using the "Thicken Wire" tool (or in Python: XCoreModeling.ThickenWire).

Since this seems to be an important question, maybe somebody who has a working script to assign homogeneous anisotropic conductivities can share it here as a reference.

[montanaro]

Something to get started to create an anisotropic conductivity tensor for S4L from solving a Poisson equation with boundary conditions (running an EM simulation and using its potential). Once the field is created in the analysis tab it can be used as a part of a new EM simulation for any material selected to have inhomogeneous and anisotropic properties.

import numpy as np
import s4l_v1.document as document
import s4l_v1.model as model
import s4l_v1.simulation.emlf as emlf
from s4l_v1 import ReleaseVersion, Unit
from s4l_v1.model import Vec3, Translation, Rotation
import s4l_v1.units as units
import s4l_v1.analysis as analysis
import XCoreModeling

def CreateMockAnisotropicConductivity(sim):
	# Transversal and longitudinal sigmas
	sigma1 = 3.0
	sigma2 = 1.0

	# here conductivity has 6 component (anisotropic tensor)
	results = sim.Results()
	p = results["Overall Field"].Outputs["EM Potential(x,y,z,f0)"]
	p.Update()
	e = results["Overall Field"].Outputs["EM E(x,y,z,f0)"]
	e.Update()

	e_field = e.Data.Field(0)
	nc = e.Data.Grid.NumberOfCells
	print "E Field shape", e_field.shape
	print "Number of Cells", nc
	# E Field shape (nc, 3) - Cell centered

	# Calculate sigma on cell centers
	print "Calculate anisotropic sigma values"
	sigma_wm = np.zeros([nc,6]) # 6 unique values
	# tensor component indexing for solver
	iu3 = (np.array([0, 1, 2, 0, 1, 0], dtype=np.int64), np.array([0, 1, 2, 1, 2, 2], dtype=np.int64))
	for idx, val in enumerate(e_field):
		tmp = np.zeros([3,3], dtype=np.complex)
		tmp[0] = val	
		if (not np.isnan(val).any() and np.inner(val, val) != 0):
			P = np.dot(tmp.T, tmp) / np.inner(val,val)
		else:
			P = np.full((3, 3), 0, dtype=np.complex)
		sigma_tensor = sigma1*P + sigma2*( np.identity(P.shape[0]) - P ) 
		sigma_wm[idx] = sigma_tensor[iu3].real
	
	sigma_swapped = sigma_wm[:,[0,1,2,3,5,4]] # Fix bug where yz and xz were swapped
	
	tensor_conductivity = analysis.core.DoubleFieldData()
	tensor_conductivity.NumberOfSnapshots = 1
	tensor_conductivity.NumberOfComponents = 6
	tensor_conductivity.Grid = e.Data.Grid.Clone()
	tensor_conductivity.ValueLocation = analysis.core.eValueLocation.kCellCenter
	tensor_conductivity.SetField(0, sigma_swapped) #sigma_wm) # Fix bug where yz and xz were swapped
	tensor_conductivity.Quantity.Name = 'Anisotropic Conductivity From E Field'
	tensor_conductivity.Quantity.Unit = analysis.core.Unit('S/m')
	
	producer_t = analysis.core.TrivialProducer()
	producer_t.SetDataObject(tensor_conductivity)
	producer_t.Name = 'Anisotropic Conductivity'
	document.AllAlgorithms.Add(producer_t)