Commit 75f98556 authored by Emmanuel Roubin's avatar Emmanuel Roubin
Browse files

[skip-ci] fix connectivity matrix for the projection

parent 982a0bb5
Pipeline #63149 skipped
# -*- coding: utf-8 -*-
"""
Project a scalar field onto an unstructured mesh
================================================
This example shows how to generate a correlated random field project it onto an unstructured mesh
"""
# sphinx_gallery_thumbnail_number = 1
import spam.excursions
import spam.mesh
import matplotlib.pyplot as plt
#####################################
# Create the random field
#####################################
# Generate a realisation of the correlated (Gaussian centered) Random Field in a cube of size 1x1x1
covariance = {'type': 'stable', 'alpha': 2.0, 'variance': 1.0, 'correlation_length': 0.5}
nNodes = 100 # discretisation of the field
realisation = spam.excursions.simulateRandomField(nNodes=nNodes, covariance=covariance, dim=3)
plt.figure()
plt.imshow(realisation[0, :, :])
plt.title("2D Slice of the field")
plt.show()
#####################################
# Generate a cylindrical mesh
#####################################
# Set the geometrical properties of the mesh
centre = [0.5, 0.5]
radius = 0.5
height = 1.
lcar = 0.1 # average distance between nodes
points, connectivity = spam.mesh.createCylinder(centre, radius, height, lcar, vtkFile="cylinder")
# points, connectivity = spam.mesh.createCuboid([1, 1, 1], lcar, vtkFile="cube")
#######################################
# Projection
#######################################
# Set the geometrical properties of the mesh
field = {
"origin": [0, 0, 0],
"lengths": [1, 1, 1],
"nCells": [99, 99, 99],
"values": [realisation]
}
thresholds = [0.2]
connectivity, materials = spam.mesh.projectField({"points": points, "cells": connectivity}, field, thresholds=thresholds, writeConnectivity='spam_projection', vtkMesh="Hello")
#######################################
# Returns the connectivity matrix and the material repartition
print(f"Connectivity matrix: {connectivity.shape} -> (node 1, node 2, node 3, node 4)")
print(f"Materials: {materials.shape} -> (phase number, sub volume, orientation x, orientation y, orientation z)")
#######################################
# Notes
#######################################
# One can generate several fields and set them in the list ``field["values"]`` line 46.
# The same number of thresholds have to be set on line 48.
......@@ -481,7 +481,7 @@ def slicePadded(im, startStop, createMask=False, padValue=0):
mask : 3D numpy array of bools
The 3D mask, only returned if createMask is True
"""
startStop = numpy.array(startStop).astype(numpy.int)
startStop = numpy.array(startStop).astype(int)
assert (startStop[1]>startStop[0]), "spam.helpers.slicePadded(): Zmax should be bigger than Zmin"
assert (startStop[3]>startStop[2]), "spam.helpers.slicePadded(): Ymax should be bigger than Ymin"
......
......@@ -439,7 +439,7 @@ def projectField(mesh, fields, thresholds=[0.0], nSkip=1, writeConnectivity=None
-------
(connectivity, materials): (nElem x 4 numpy array, nElem x 5 numpy array)
**connectivity**: the classical connectivity matrix
**materials**: for each element ``elemType, interX, interY, interZ, subVol`` (see outputFile for details).
**materials**: for each element ``elemType, subVol, interX, interY, interZ`` (see outputFile for details).
"""
# init projmorpho
......@@ -456,8 +456,43 @@ def projectField(mesh, fields, thresholds=[0.0], nSkip=1, writeConnectivity=None
points = mesh['points']
cells = mesh['cells']
# DEBUG: Artificially add unused node
# points = numpy.insert(points, 10, [-1.0,-1.0,-1.0], axis=0)
# for i, node_number_x4 in enumerate(cells):
# for j, node_number in enumerate(node_number_x4):
# if cells[i, j] >= 10:
# cells[i, j] += 1
# test if unused nodes
cells_set = set(cells.ravel()) # ordered numbering in connectivity
n_points_used = len(cells_set)
if n_points_used != points.shape[0]:
print("Deleting points before projection")
# removing unused nodes
points_to_delete = []
n_deleted = 0
for new, old in enumerate(cells_set):
# check if holes in the connectivity matrix
if new != old - n_deleted:
points_to_delete.append(new + n_deleted)
n_deleted += 1
# print(new, old, n_deleted, old - n_deleted)
points = numpy.delete(points, points_to_delete, axis=0)
print("Renumbering connectivity matrix")
# renumbering connectivity matrix accordingly
new_node_numbering = {old: new + 1 for new, old in enumerate(cells_set)}
for i, node_number_x4 in enumerate(cells):
for j, node_number in enumerate(node_number_x4):
if new_node_numbering[cells[i, j]] != cells[i, j]:
cells[i, j] = new_node_numbering[cells[i, j]]
del cells_set
# check if cell number start by 1 or 0
if cells.min() == 0:
print("Shift connectivity by 1 to start at 1")
cells += 1
pr.setMesh(points, cells.ravel())
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment