<commit_before>"""
Generates 3-D Left and Right ventricles mesh starting from modified sphere shell mesh.
"""

import math
from mapclientplugins.meshgeneratorstep.meshtypes.meshtype_3d_sphereshell1 import MeshType_3d_sphereshell1
from mapclientplugins.meshgeneratorstep.utils.eft_utils import *
from mapclientplugins.meshgeneratorstep.utils.zinc_utils import *
from mapclientplugins.meshgeneratorstep.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
from opencmiss.zinc.element import Element, Elementbasis, Elementfieldtemplate
from opencmiss.zinc.field import Field
from opencmiss.zinc.node import Node

class MeshType_3d_heartventricles1:
    '''
    classdocs
    '''
    @staticmethod
    def getName():
        return '3D Heart Ventricles 1'

    @staticmethod
    def getDefaultOptions():
        return {
            'Number of elements up' : 4,
            'Number of elements around' : 10,
            'Number of elements across septum' : 4,
            'Number of elements below septum' : 2,
            'Number of elements through LV wall' : 1,
            'LV wall thickness' : 0.15,
            'LV wall thickness ratio apex' : 0.5,
            'RV free wall thickness' : 0.05,
            'RV width' : 0.2,
            'Length ratio' : 2.0,
            'Element length ratio equator/apex' : 1.0,
            'Use cross derivatives' : False
        }

    @staticmethod
    def getOrderedOptionNames():
        return [
            'Number of elements up',
            'Number of elements around',
            'Number of elements across septum',
            'Number of elements below septum',
            'Number of elements through LV wall',
            'LV wall thickness',
            'LV wall thickness ratio apex',
            'RV free wall thickness',
            'RV width',
            'Length ratio',
            'Element length ratio equator/apex'
        ]

    @staticmethod
    def checkOptions(options):
        if options['Number of elements up'] < 4:
            options['Number of elements up'] = 4
        if options['Number of elements around'] < 4:
            options['Number of elements around'] = 4
        if options['Number of elements across septum'] < 2:
            options['Number of elements across septum'] = 2
        elif options['Number of elements across septum'] > (options['Number of elements around'] - 2):
            options['Number of elements across septum'] = options['Number of elements around'] - 2
        if options['Number of elements below septum'] < 2:
            options['Number of elements below septum'] = 2
        elif options['Number of elements below septum'] > (options['Number of elements up'] - 1):
            options['Number of elements below septum'] = options['Number of elements up'] - 1
        if (options['Number of elements through LV wall'] < 1) :
            options['Number of elements through LV wall'] = 1
        if options['LV wall thickness'] < 0.0:
            options['LV wall thickness'] = 0.0
        elif options['LV wall thickness'] > 0.5:
            options['LV wall thickness'] = 0.5
        if options['LV wall thickness ratio apex'] < 0.0:
            options['LV wall thickness ratio apex'] = 0.0
        if options['RV free wall thickness'] < 0.0:
            options['RV free wall thickness'] = 0.0
        if options['RV width'] < 0.0:
            options['RV width'] = 0.0
        if options['Length ratio'] < 1.0E-6:
            options['Length ratio'] = 1.0E-6
        if options['Element length ratio equator/apex'] < 1.0E-6:
            options['Element length ratio equator/apex'] = 1.0E-6

    @staticmethod
    def generateMesh(region, options):
        """
        :param region: Zinc region to define model in. Must be empty.
        :param options: Dict containing options. See getDefaultOptions().
        :return: None
        """
        elementsCountUp = options['Number of elements up']
        elementsCountAround = options['Number of elements around']
        elementsCountAcrossSeptum = options['Number of elements across septum']
        elementsCountBelowSeptum = options['Number of elements below septum']
        elementsCountThroughLVWall = options['Number of elements through LV wall']
        LVWallThickness = options['LV wall thickness']
        LVWallThicknessRatioApex = options['LV wall thickness ratio apex']
        RVWidth = options['RV width']
        RVFreeWallThickness = options['RV free wall thickness']
        useCrossDerivatives = options['Use cross derivatives']

        # generate a half sphere shell which will be edited
        sphereShellOptions = MeshType_3d_sphereshell1.getDefaultOptions()
        sphereShellOptions['Number of elements up'] = elementsCountUp*2
        sphereShellOptions['Number of elements around'] = elementsCountAround
        sphereShellOptions['Number of elements through wall'] = elementsCountThroughLVWall
        sphereShellOptions['Exclude top rows'] = elementsCountUp
        sphereShellOptions['Wall thickness'] = LVWallThickness
        sphereShellOptions['Wall thickness ratio apex'] = LVWallThicknessRatioApex
        sphereShellOptions['Length ratio'] = options['Length ratio']
        sphereShellOptions['Element length ratio equator/apex'] = options['Element length ratio equator/apex']
        MeshType_3d_sphereshell1.generateMesh(region, sphereShellOptions)

        fm = region.getFieldmodule()
        fm.beginChange()
        # find the coordinates field created for the sphere shell
        coordinates = fm.findFieldByName('coordinates').castFiniteElement()

        nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        nodetemplate = nodes.createNodetemplate()
        nodetemplate.defineField(coordinates)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)

        mesh = fm.findMeshByDimension(3)

        tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives)
        tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE)

        eft = tricubichermite.createEftBasic()
        elementtemplate = mesh.createElementtemplate()
        elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE)
        elementtemplate.defineField(coordinates, -1, eft)

        cache = fm.createFieldcache()

        edgeOffsetInner = RVWidth*0.125
        edgeOffsetOuter = edgeOffsetInner*2.0 + RVFreeWallThickness
        crossAngle = math.pi/12
        sinCrossAngle = math.sin(crossAngle)
        cosCrossAngle = math.cos(crossAngle)
        edgeRotateCrossAngle = math.pi/24
        sinEdgeRotateCrossAngle = math.sin(edgeRotateCrossAngle)
        cosEdgeRotateCrossAngle = math.cos(edgeRotateCrossAngle)

        # create RV nodes and modify adjoining LV nodes
        elementsCountUpRV = elementsCountUp - elementsCountBelowSeptum
        now = 1 + elementsCountUp*elementsCountAround
        nodeIdentifier = 1 + 2*elementsCountThroughLVWall*now
        rv_nidsInner = []
        rv_nidsOuter = []
        for n3 in range(1, -1, -1):  # only 1 element through RV free wall

            onInside = n3 == 0

            for n2 in range(elementsCountUpRV + 1):

                onBottomEdge = (n2 == 0)

                for n1 in range(elementsCountAcrossSeptum + 1):
<commit_msg>Add left ventricle wall thinning at base<commit_after>"""
Generates 3-D Left and Right ventricles mesh starting from modified sphere shell mesh.
"""

import math
from mapclientplugins.meshgeneratorstep.meshtypes.meshtype_3d_sphereshell1 import MeshType_3d_sphereshell1
from mapclientplugins.meshgeneratorstep.utils.eft_utils import *
from mapclientplugins.meshgeneratorstep.utils.zinc_utils import *
from mapclientplugins.meshgeneratorstep.utils.eftfactory_tricubichermite import eftfactory_tricubichermite
from opencmiss.zinc.element import Element, Elementbasis, Elementfieldtemplate
from opencmiss.zinc.field import Field
from opencmiss.zinc.node import Node

class MeshType_3d_heartventricles1:
    '''
    classdocs
    '''
    @staticmethod
    def getName():
        return '3D Heart Ventricles 1'

    @staticmethod
    def getDefaultOptions():
        return {
            'Number of elements up' : 4,
            'Number of elements around' : 10,
            'Number of elements across septum' : 4,
            'Number of elements below septum' : 2,
            'Number of elements through LV wall' : 1,
            'LV wall thickness' : 0.15,
            'LV wall thickness ratio apex' : 0.5,
            'LV wall thickness ratio base': 0.5,
            'RV free wall thickness' : 0.05,
            'RV width' : 0.2,
            'Length ratio' : 2.0,
            'Element length ratio equator/apex' : 1.0,
            'Use cross derivatives' : False
        }

    @staticmethod
    def getOrderedOptionNames():
        return [
            'Number of elements up',
            'Number of elements around',
            'Number of elements across septum',
            'Number of elements below septum',
            'Number of elements through LV wall',
            'LV wall thickness',
            'LV wall thickness ratio apex',
            'LV wall thickness ratio base',
            'RV free wall thickness',
            'RV width',
            'Length ratio',
            'Element length ratio equator/apex'
        ]

    @staticmethod
    def checkOptions(options):
        if options['Number of elements up'] < 4:
            options['Number of elements up'] = 4
        if options['Number of elements around'] < 4:
            options['Number of elements around'] = 4
        if options['Number of elements across septum'] < 2:
            options['Number of elements across septum'] = 2
        elif options['Number of elements across septum'] > (options['Number of elements around'] - 2):
            options['Number of elements across septum'] = options['Number of elements around'] - 2
        if options['Number of elements below septum'] < 2:
            options['Number of elements below septum'] = 2
        elif options['Number of elements below septum'] > (options['Number of elements up'] - 1):
            options['Number of elements below septum'] = options['Number of elements up'] - 1
        if (options['Number of elements through LV wall'] < 1) :
            options['Number of elements through LV wall'] = 1
        if options['LV wall thickness'] < 0.0:
            options['LV wall thickness'] = 0.0
        elif options['LV wall thickness'] > 0.5:
            options['LV wall thickness'] = 0.5
        if options['LV wall thickness ratio apex'] < 0.0:
            options['LV wall thickness ratio apex'] = 0.0
        if options['LV wall thickness ratio base'] < 0.0:
            options['LV wall thickness ratio base'] = 0.0
        if options['RV free wall thickness'] < 0.0:
            options['RV free wall thickness'] = 0.0
        if options['RV width'] < 0.0:
            options['RV width'] = 0.0
        if options['Length ratio'] < 1.0E-6:
            options['Length ratio'] = 1.0E-6
        if options['Element length ratio equator/apex'] < 1.0E-6:
            options['Element length ratio equator/apex'] = 1.0E-6

    @staticmethod
    def generateMesh(region, options):
        """
        :param region: Zinc region to define model in. Must be empty.
        :param options: Dict containing options. See getDefaultOptions().
        :return: None
        """
        elementsCountUp = options['Number of elements up']
        elementsCountAround = options['Number of elements around']
        elementsCountAcrossSeptum = options['Number of elements across septum']
        elementsCountBelowSeptum = options['Number of elements below septum']
        elementsCountThroughLVWall = options['Number of elements through LV wall']
        LVWallThickness = options['LV wall thickness']
        LVWallThicknessRatioApex = options['LV wall thickness ratio apex']
        LVWallThicknessRatioBase = options['LV wall thickness ratio base']
        RVWidthTop = options['RV width']
        RVFreeWallThickness = options['RV free wall thickness']
        useCrossDerivatives = options['Use cross derivatives']

        # generate a half sphere shell which will be edited
        sphereShellOptions = MeshType_3d_sphereshell1.getDefaultOptions()
        sphereShellOptions['Number of elements up'] = elementsCountUp*2
        sphereShellOptions['Number of elements around'] = elementsCountAround
        sphereShellOptions['Number of elements through wall'] = elementsCountThroughLVWall
        sphereShellOptions['Exclude top rows'] = elementsCountUp
        sphereShellOptions['Wall thickness'] = LVWallThickness
        sphereShellOptions['Wall thickness ratio apex'] = LVWallThicknessRatioApex
        sphereShellOptions['Length ratio'] = options['Length ratio']
        sphereShellOptions['Element length ratio equator/apex'] = options['Element length ratio equator/apex']
        MeshType_3d_sphereshell1.generateMesh(region, sphereShellOptions)

        fm = region.getFieldmodule()
        fm.beginChange()
        # find the coordinates field created for the sphere shell
        coordinates = fm.findFieldByName('coordinates').castFiniteElement()

        nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        nodetemplate = nodes.createNodetemplate()
        nodetemplate.defineField(coordinates)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)

        mesh = fm.findMeshByDimension(3)

        cache = fm.createFieldcache()

        if LVWallThicknessRatioBase != 1.0:
            # make LV walls thinner at base
            # get inside node at middle of RV
            now = 1 + elementsCountUp*elementsCountAround
            midRVnid = now - elementsCountAround + 2 + (elementsCountAcrossSeptum // 2)
            #print('midRVnid', midRVnid)
            midRVnode = nodes.findNodeByIdentifier(midRVnid)
            cp_coordinates = fm.createFieldCoordinateTransformation(coordinates)
            cp_coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_CYLINDRICAL_POLAR)
            radius = fm.createFieldComponent(cp_coordinates, 1)
            theta = fm.createFieldComponent(cp_coordinates, 2)
            z = fm.createFieldComponent(cp_coordinates, 3)
            cache.setNode(midRVnode)
            #result, cp = cp_coordinates.evaluateReal(cache, 3)
            #coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3)
            result, innerRadius = radius.evaluateReal(cache, 1)
            #print('innerRadius', innerRadius)
            ir = fm.createFieldConstant([innerRadius*0.9999])
            radius_gt_ir = fm.createFieldGreaterThan(radius, ir)
            radius_minus_ir = fm.createFieldSubtract(radius, ir)
            thickness_scale = fm.createFieldConstant([LVWallThicknessRatioBase])
            delta_radius = fm.createFieldMultiply(radius_minus_ir, thickness_scale)
            new_radius = fm.createFieldAdd(ir, delta_radius)
            new_cp_coordinates = fm.createFieldConcatenate([new_radius, theta, z])
            new_cp_coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_CYLINDRICAL_POLAR)
            new_coordinates = fm.createFieldCoordinateTransformation(new_cp_coordinates)
            new_coordinates.setCoordinateSystemType(Field.COORDINATE_SYSTEM_TYPE_RECTANGULAR_CARTESIAN)

            baseNodeGroupField = fm.createFieldNodeGroup(nodes)
            just_under_zero = fm.createFieldConstant([-0.0001])
            isBase = fm.createFieldGreaterThan(z, just_under_zero)
            baseNodesetGroup = baseNodeGroupField.getNodesetGroup()
            baseNodesetGroup.addNodesConditional(isBase)
            #print('baseNodesetGroup.getSize()', baseNodesetGroup.getSize())
            fieldassignment = coordinates.createFieldassignment(new_coordinates)
            result = fieldassignment.setNodeset(baseNodesetGroup)
            #print('fieldassignment.setNodeset', result)
            fieldassignment.assign()

        tricubichermite = eftfactory_tricubichermite(mesh, useCrossDerivatives)
        tricubicHermiteBasis = fm.createElementbasis(3, Elementbasis.FUNCTION_TYPE_CUBIC_HERMITE)

        eft = tricubichermite.createEftBasic()
        elementtemplate = mesh.createElementtemplate()
        elementtemplate.setElementShapeType(Element.SHAPE_TYPE_CUBE)
        elementtemplate.defineField(coordinates, -1, eft)

        crossAngle = math.pi/8
        sinCrossAngle = math.sin(crossAngle)
        cosCrossAngle = math.cos(crossAngle)
        edgeRotateCrossAngle = math.pi/16
        sinEdgeRotateCrossAngle = math.sin(edgeRotateCrossAngle)
        cosEdgeRotateCrossAngle = math.cos(edgeRotateCrossAngle)

        # create RV nodes and modify adjoining LV nodes
        elementsCountUpRV = elementsCountUp - elementsCountBelowSeptum
        now = 1 + elementsCountUp*elementsCountAround
        nodeIdentifier = 1 + 2*elementsCountThroughLVWall*now
        rv_nidsInner = []
        rv_nidsOuter = []
        baseExtraRVWidth = LVWallThickness*(1.0 - LVWallThicknessRatioBase)
        for n3 in range(1, -1, -1):  # only 1 element through RV free wall

            onInside = n3 == 0

            for n2 in range(elementsCountUpRV + 1):

                #if n2 > (elementsCountUpRV - elementsCountUpBase):
                #    theta = math.pi*0.5
                #else:
                #    theta = (n2 / (elementsCountUpRV - elementsCountUpBase))*math.pi*0.5
                #theta = (n2 / elementsCountUpRV)*math.pi*0.5
                RVWidth = RVWidthTop  #*math.sin(theta)
                #RVWidth = RVWidthTop

                edgeOffsetInner = RVWidth*0.125
                if n2 == elementsCountUpRV:
                    RVWidth += baseExtraRVWidth
                    edgeOffsetInner = edgeOffsetInner + baseExtraRVWidth*0.5

                onBottomEdge = (n2 == 0)

                for n1 in range(elementsCountAcrossSeptum + 1):

{'fraction_non_alphanumeric': 0.08455346047868605, 'fraction_numerical': 0.015127775676697822, 'mean_word_length': 2.725845410628019, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2018049', 'n_tokens_mistral': 5329, 'n_tokens_neox': 4950, 'n_words': 1337}