asdeasoft.net

Code: Select all

from PyMpc import *
from PyMpc import MpcOdbVirtualResult as vr
import math

App.clearTerminal()
doc = App.postDocument()
pg = doc.activePlotGroup

# clear current selection
doc.scene.unselectAll()

# get a result and evaluate it to get a mesh
database_id = 1
db = doc.getDatabase(database_id)
U = db.getNodalResult('Displacement')
all_stages = db.getStageIDs()
last_stage = all_stages[-1]
all_steps = db.getStepIDs(last_stage)
last_step = all_steps[-1]
opt = MpcOdbVirtualResultEvaluationOptions()
opt.stage = last_stage
opt.step = last_step
field = U.evaluate(opt)
mesh = field.mesh

# select all elements according to the position of their centroid.
for ele_id, ele in mesh.elements.items():
	c = ele.computeCenter()
	# here you can implement your logic.
	# In this example I want all elements at a positive X
	if c.x > 0:
		# select this element in all plots of the current plot group
		for plot_id, plot in pg.plots.items():
			# get the selection data for this plot, or create it
			selection_data = doc.scene.plotSelection.get(plot_id, None)
			if selection_data is None:
				selection_data = MpcOdpSelectionData()
				doc.scene.plotSelection[plot_id] = selection_data
			# add this element
			selection_data.info.elements.insert(ele_id)
# done
doc.scene.updateSelectionVisualRepresentationAndGraphics()
App.updateActiveView()

Code: Select all

from PyMpc import *

# USER SETTINGS
DATABASE_ID = 1
ARROW_SCALE = 1.0

App.clearTerminal()
doc = App.postDocument()
doc.clearCustomDrawableEntities()

# get active plot group
group = doc.activePlotGroup

# allocate a shape for plotting
vrep = FxShape()
vrep_edge = None

# process all plots in this plot group
for _, plot in group.plots.items():
	# get source DB, stage and step indices
	db = doc.getDatabase(plot.plotData.db)
	stage = plot.plotData.stage
	step = plot.plotData.step
	# get result
	P = db.getNodalResult('Pressure')
	# evaluate it
	opt = MpcOdbVirtualResultEvaluationOptions()
	opt.stage = stage
	opt.step = step
	PField = P.evaluate(opt)
	# now we have a scalar field of the water pore pressure.
	# we can compute its spatial derivative DP = dP/dX (a vector field)
	# -DP will give the path
	N = Math.vec() # pre-allocate a vector for shape functions
	dN = Math.mat() # pre-allocate a matrix for shape function local derivatives
	J = Math.mat() # pre-allocate a matrix for jacobian
	for _, ele in PField.mesh.elements.items():
		# compute the location of the center in natural coords
		p0 = MpcIntegrationPoint(0.0, 0.0)
		for ip in ele.integrationRule.integrationPoints:
			p0.x += ip.x
			p0.y += ip.y
			p0.w += ip.w
		p0.x /= ele.numberOfIntegrationPoints()
		p0.y /= ele.numberOfIntegrationPoints()
		# evaluate local derivatives and jacobian at center
		ele.jacobianAt(p0, dN, J)
		# handle the 2D case
		if J.cols == 2:
			dN0 = dN.col(0)
			dN1 = dN.col(1)
			dN = Math.mat(dN.rows, 3, 0.0)
			for i in range(dN.rows):
				dN[i, 0] = dN0[i]
				dN[i, 1] = dN1[i]
			J0 = J.col(0)
			J1 = J.col(1)
			J = Math.mat(3,3, 0.0)
			for i in range(3):
				J[i, 0] = J0[i]
				J[i, 1] = J1[i]
			J[2, 2] = 1.0
		# compute cartesian derivatives
		dNdX = dN * J.inverse()
		# collect pressure values at the nodes
		PNodal = Math.vec([ PField[MpcOdbResultField.node(node.id), 0] for node in ele.nodes ])
		# compute the gradient of the scalar Pressure at the center
		dPdX = Math.vec3(*(dNdX.transpose() * PNodal))
		dPdX *= -1.0 # negate it, we want the direction of decreasing pressure
		# update VREP
		center = ele.computeCenter()
		length = dPdX.norm() * ARROW_SCALE
		axis = dPdX.normalized()
		radius_ratio = 0.04
		num_seg = 6
		double_arrow = False
		vrep_arrow = FxShapeFactory.arrowLine(center, axis, length, radius_ratio, num_seg, double_arrow)
		if vrep_edge is None:
			vrep_edge = vrep_arrow
		else:
			vrep_edge.merge(vrep_arrow)

# add the VREP to the scene
vrep.edges.append(vrep_edge)
vrep_material = FxMaterial()
vrep_material.edgeColor = FxColor(.6,0,0) # R,G,B in range[0,1]
vrep_material.lineWidth = 2.0
vrep.material = vrep_material
doc.addCustomDrawableEntity(vrep)
App.updateActiveView()