diff --git a/Assembler.py b/Assembler.py index 65f3cad..621371f 100755 --- a/Assembler.py +++ b/Assembler.py @@ -9,8 +9,8 @@ ## Copyright (c) Simone Manini, Luca Antiga. All rights reserved. ## See LICENCE file for details. -## This software is distributed WITHOUT ANY WARRANTY; without even -## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR ## PURPOSE. See the above copyright notices for more information. ## Developed with support from the EC FP7/2007-2013: ARCH, Project n. 224390 @@ -29,7 +29,7 @@ class Assembler(object): Assemble: a method for building boundary conditions vectors and assembling each local Zero, First and Second Order matrix into global Zero, First and Second Order matrix. ''' - + def __init__(self): ''' Constructor @@ -48,13 +48,13 @@ def __init__(self): self.SecondOrderGlobalMatrix = None self.Evaluator = None self.Initialized = False - + def SetNetworkMesh(self, networkMesh): ''' Setting NetworkMesh ''' self.NetworkMesh = networkMesh - + def SetBoundaryConditions(self, boundaryConditions): ''' Setting BoundaryConditions @@ -66,7 +66,7 @@ def GetNumberOfGlobalDofs(self): This method returns number of global degrees of freedom ''' return self.DofMap.NumberOfGlobalDofs - + def AssembleBoundaryConditions(self, simulationContext): ''' This method assembles arrays of prescribed pressures @@ -74,33 +74,35 @@ def AssembleBoundaryConditions(self, simulationContext): self.DofMap = DofMap() self.DofMap.SetNetworkMesh(self.NetworkMesh) self.DofMap.Build() - + # Searching for prescribed pressure output/input numberOfElements = 0 - for element in self.NetworkMesh.Elements: + for element in self.NetworkMesh.Elements: for dof in element.GetExternalPressureLocalDofs(): numberOfElements+=1 if self.BoundaryConditions.OutP is not None: - numberOfElements+=1 + numberOfElements+=len(self.BoundaryConditions.PressureOut) if self.BoundaryConditions.InP is not None: numberOfElements+=1 - - + # Setting Transmural Pressures for windkessel elements and wave propagation elements PrescribedPressures = zeros((numberOfElements,2)) done = 0 i = 0 for element in self.NetworkMesh.Elements: - for dof in element.GetExternalPressureLocalDofs(): - if self.BoundaryConditions.OutP is not None: - if element.Id == self.BoundaryConditions.elementOut.Id: - if done == 0: - value1 = self.DofMap.DofMap[(self.BoundaryConditions.elementOut.Id,self.BoundaryConditions.elementOut.GetLocalDof(int(self.BoundaryConditions.NodeOut)))] - value2 = self.BoundaryConditions.OutP - PrescribedPressures[i,0] = value1 - PrescribedPressures[i,1] = value2 - done = 1 - i+=1 + for dof in element.GetExternalPressureLocalDofs(): + + for el, elProps in self.BoundaryConditions.PressureOut.iteritems(): + if element.Id == el.Id: + if done < len(self.BoundaryConditions.PressureOut): + value1 = self.DofMap.DofMap[(el.Id,el.GetLocalDof(int(elProps['node'])))] + value2 = elProps['value'] + if value1 not in PrescribedPressures: + PrescribedPressures[i,0] = value1 + PrescribedPressures[i,1] = value2 + done+=1 + i+=1 + if self.BoundaryConditions.InP is not None: if element.Id == self.BoundaryConditions.elementIn.Id: if done == 0: @@ -110,13 +112,13 @@ def AssembleBoundaryConditions(self, simulationContext): PrescribedPressures[i,1] = value2 done = 1 i+=1 - + value1 = self.DofMap.DofMap[(element.Id,dof)] value2 = self.BoundaryConditions.PressureValues[element.Id] PrescribedPressures[i,0] = value1 - PrescribedPressures[i,1] = value2 + PrescribedPressures[i,1] = value2 i+=1 - + self.BoundaryConditions.SetSimulationContext(simulationContext) for el in self.BoundaryConditions.elementFlow: self.FlowDof[el.Id] = self.DofMap.DofMap[(el.Id,el.GetLocalDof(int(self.BoundaryConditions.NodeFlow[el.Id])))] @@ -129,7 +131,7 @@ def AssembleInit(self, simulationContext, evaluator): Pressures GlobalDofs and Input Flow Vector from BoundaryConditions (Prescribed Pressures matrix, each row ---> Dof, Value). Zero, First and Second Order Global Matrix from Local Matrix of each element. ''' - + #Assembling global matrices from local matrices. self.LinearZeroOrderGlobalMatrix = zeros((self.DofMap.NumberOfGlobalDofs, self.DofMap.NumberOfGlobalDofs)) self.LinearFirstOrderGlobalMatrix = zeros((self.DofMap.NumberOfGlobalDofs, self.DofMap.NumberOfGlobalDofs)) @@ -137,13 +139,13 @@ def AssembleInit(self, simulationContext, evaluator): self.ZeroOrderGlobalMatrix = zeros((self.DofMap.NumberOfGlobalDofs, self.DofMap.NumberOfGlobalDofs)) self.FirstOrderGlobalMatrix = zeros((self.DofMap.NumberOfGlobalDofs, self.DofMap.NumberOfGlobalDofs)) self.SecondOrderGlobalMatrix = zeros((self.DofMap.NumberOfGlobalDofs, self.DofMap.NumberOfGlobalDofs)) - + nonLinear = False - #Building Global Linear Matrices + #Building Global Linear Matrices for element in self.DofMap.NetworkMesh.Elements: - if element.Initialized == False: + if element.Initialized == False: element.Initialize(simulationContext) - if element.nonLinear == False: + if element.nonLinear == False: element.InputParameters(evaluator) zeroOrderMatrix = element.GetZeroOrderMatrix() firstOrderMatrix = element.GetFirstOrderMatrix() @@ -152,34 +154,34 @@ def AssembleInit(self, simulationContext, evaluator): rowGlobalDof = self.DofMap.GetDof(element.Id,rowLocalDof) for columnLocalDof in element.dof: columnGlobalDof = self.DofMap.GetDof(element.Id,columnLocalDof) - self.LinearZeroOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += zeroOrderMatrix[rowLocalDof,columnLocalDof] + self.LinearZeroOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += zeroOrderMatrix[rowLocalDof,columnLocalDof] self.LinearFirstOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += firstOrderMatrix[rowLocalDof,columnLocalDof] self.LinearSecondOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += secondOrderMatrix[rowLocalDof,columnLocalDof] else: nonLinear = True - - if nonLinear == True: + + if nonLinear == True: self.ZeroOrderGlobalMatrix,self.FirstOrderGlobalMatrix,self.SecondOrderGlobalMatrix = \ self.Assemble(simulationContext, evaluator, self.LinearZeroOrderGlobalMatrix, self.LinearFirstOrderGlobalMatrix, self.LinearSecondOrderGlobalMatrix) - + else: self.ZeroOrderGlobalMatrix[:,:] = self.LinearZeroOrderGlobalMatrix[:,:] self.FirstOrderGlobalMatrix[:,:] = self.LinearFirstOrderGlobalMatrix[:,:] self.SecondOrderGlobalMatrix[:,:] = self.LinearSecondOrderGlobalMatrix[:,:] - + return self.LinearZeroOrderGlobalMatrix, self.LinearFirstOrderGlobalMatrix, self.LinearSecondOrderGlobalMatrix - - - def Assemble(self, simulationContext, evaluator, linearZeroOrder, linearFirstOrder, linearSecondOrder): + + + def Assemble(self, simulationContext, evaluator, linearZeroOrder, linearFirstOrder, linearSecondOrder): ''' This method builds non-linear elements of the global matrices from linear matrices. ''' - + #Building non linear matrices - self.ZeroOrderGlobalMatrix[:,:] = linearZeroOrder[:,:] + self.ZeroOrderGlobalMatrix[:,:] = linearZeroOrder[:,:] self.FirstOrderGlobalMatrix[:,:] = linearFirstOrder[:,:] - self.SecondOrderGlobalMatrix[:,:] = linearSecondOrder[:,:] - + self.SecondOrderGlobalMatrix[:,:] = linearSecondOrder[:,:] + for element in self.DofMap.NetworkMesh.Elements: if element.nonLinear == True: element.Initialize(simulationContext) @@ -191,7 +193,7 @@ def Assemble(self, simulationContext, evaluator, linearZeroOrder, linearFirstOrd rowGlobalDof = self.DofMap.GetDof(element.Id,rowLocalDof) for columnLocalDof in element.dof: columnGlobalDof = self.DofMap.GetDof(element.Id,columnLocalDof) - self.ZeroOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += zeroOrderMatrix[rowLocalDof,columnLocalDof] + self.ZeroOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += zeroOrderMatrix[rowLocalDof,columnLocalDof] self.FirstOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += firstOrderMatrix[rowLocalDof,columnLocalDof] self.SecondOrderGlobalMatrix[rowGlobalDof,columnGlobalDof] += secondOrderMatrix[rowLocalDof,columnLocalDof] - return self.ZeroOrderGlobalMatrix, self.FirstOrderGlobalMatrix, self.SecondOrderGlobalMatrix \ No newline at end of file + return self.ZeroOrderGlobalMatrix, self.FirstOrderGlobalMatrix, self.SecondOrderGlobalMatrix diff --git a/BoundaryConditions.py b/BoundaryConditions.py index e735113..1867b26 100755 --- a/BoundaryConditions.py +++ b/BoundaryConditions.py @@ -9,8 +9,8 @@ ## Copyright (c) Simone Manini, Luca Antiga. All rights reserved. ## See LICENCE file for details. -## This software is distributed WITHOUT ANY WARRANTY; without even -## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR ## PURPOSE. See the above copyright notices for more information. ## Developed with support from the EC FP7/2007-2013: ARCH, Project n. 224390 @@ -36,27 +36,28 @@ class BoundaryConditions(object): GetPressure: a method for calculating transmural pressures for a specific time value. Timestep and period from SimulationContext are necessary. ReadFromXML: a method for reading Boundary Conditions XML File and - setting boundary conditions' parameters. - - ''' - + setting boundary conditions' parameters. + + ''' + def __init__(self): ''' Constructor ''' self.Id = None - self.BC = {} #PressureValue:[elements] - self.TimePressure = {} #time:PressureValue + self.BC = {} # PressureValue:[elements] + self.TimePressure = {} # time:PressureValue self.NetworkMesh = None self.SimulationContext = None self.Flow = None self.elementFlow = [] - self.NodeFlow = {} #element:NodeFlow - self.elementOut = None - self.NodeOut = None + self.NodeFlow = {} # element:NodeFlow + self.PressureOut = {} # element: {node:node, value: value} + self.elementsOut = [] # multiple elements for boundary condition pressure + self.NodesOut = [] # multiple nodes for boundary condition pressure self.elementIn = None self.NodeIn = None - self.OutP = None + self.OutsP = [] # multiple values for boundary condition pressure self.InP = None self.PressureValues = {} # dictionary of external pressures (element:value) self.InFlows = {} # dictionary of inlet flows (element:{name:name, A0:value, f_coeff:value, signal:value) @@ -68,19 +69,19 @@ def SetSimulationContext(self,simulationContext): ''' Setting SimulationContext ''' - self.SimulationContext = simulationContext - + self.SimulationContext = simulationContext + def SetNetworkMesh(self,networkMesh): ''' Setting NetworkMesh ''' - self.NetworkMesh = networkMesh - + self.NetworkMesh = networkMesh + def SetSpecificCardiacOutput(self): ''' Adapting cardiac inflow according to specific mean cardiac output value. ''' - + for data in self.InFlows.itervalues(): if data['name'] == 'heart': try: @@ -94,13 +95,13 @@ def SetSpecificCardiacOutput(self): shift = 9.18388663e-06 k =((A1+shift)/(A0+shift)) A0 = A1 - data['f_coeff'] = f_coeff*k + data['f_coeff'] = f_coeff*k + - def GetSteadyFlow(self, el, timestep, time): ''' Calculating flow as steady (mean A0 value) - ''' + ''' A0 = self.InFlows[el]['A0'] if time < (10*timestep): Flow = A0*((time/timestep)/10) @@ -108,7 +109,7 @@ def GetSteadyFlow(self, el, timestep, time): Flow = A0 self.Flow = Flow return Flow - + def GetFlow(self): ''' Calculating inlet flow (coefficients of the FFT x(t)=A0+sum(2*Ck*exp(j*k*2*pi*f*t))) @@ -124,17 +125,17 @@ def GetFlow(self): except KeyError: print "Error, Please set period in Simulation Context XML File" raise - + t = arange(0.0,period+timestep,timestep).reshape((1,ceil(period/timestep+1.0))) Cc = self.f_coeff*1.0/2.0*1e-6 Flow = zeros((1, ceil(period/timestep+1.0))) for freq in arange(0,ceil(period/timestep+1.0)): Flow[0, freq] = self.A0_v for k in arange(0,self.f_coeff.shape[0]): - Flow[0, freq] = Flow[0, freq]+real(2.0*complex(Cc[k,0],Cc[k,1])*exp(1j*(k+1)*2.0*pi*t[0,freq]/period)) - self.Flow = Flow + Flow[0, freq] = Flow[0, freq]+real(2.0*complex(Cc[k,0],Cc[k,1])*exp(1j*(k+1)*2.0*pi*t[0,freq]/period)) + self.Flow = Flow return Flow - + def GetTimeFlow(self, el, time): ''' Calculating inlet flow (coefficients of the FFT x(t)=A0+sum(2*Ck*exp(j*k*2*pi*f*t))) @@ -146,7 +147,7 @@ def GetTimeFlow(self, el, time): except KeyError: print "Error, Please set period in Simulation Context XML File" raise - + try: signal = self.InFlows[el]['signal'] try: @@ -168,7 +169,7 @@ def GetTimeFlow(self, el, time): Flow += real(2.0*complex(Cc[k,0],Cc[k,1])*exp(1j*(k+1)*2.0*pi*time/period)) self.Flow = Flow return Flow - + def GetPressure(self,time, entity = None): ''' Calculating transmural pressures for a specific time value. @@ -176,10 +177,10 @@ def GetPressure(self,time, entity = None): TimedPressures = {} if entity is None: time = str(time) - for mesh, timepress in self.PressureValues.iteritems(): + for mesh, timepress in self.PressureValues.iteritems(): try: if timepress.has_key(time): - TimedPressures[mesh] = timepress[time] + TimedPressures[mesh] = timepress[time] except AttributeError: TimedPressures[mesh] = timepress if entity is not None: @@ -191,11 +192,11 @@ def GetPressure(self,time, entity = None): if el.Id == mesh: try: if timepress.has_key(time): - TimedPressures[mesh] = timepress[time] + TimedPressures[mesh] = timepress[time] except AttributeError: - TimedPressures[mesh] = timepress + TimedPressures[mesh] = timepress return TimedPressures - + def ReadFromXML(self, xmlBcpath, xsdBcpath=None): ''' This method reads Boundary Conditions XML File. @@ -209,7 +210,7 @@ def ReadFromXML(self, xmlBcpath, xsdBcpath=None): LXMLError() lxml = False from xml.etree import ElementTree as etree - + if lxml: if not xsdBcpath: NoXSDWarning() @@ -227,7 +228,7 @@ def ReadFromXML(self, xmlBcpath, xsdBcpath=None): xmlschema.assert_(docbc) print "Boundary Conditions Xml File has been validated." break - except AssertionError: + except AssertionError: XMLValidationError(xmlschema) docbcfile = open(xmlBcpath) @@ -237,7 +238,7 @@ def ReadFromXML(self, xmlBcpath, xsdBcpath=None): self.Id = bcgraph_dict['id'] if self.Id != self.NetworkMesh.Id: raise XMLIdError() - + for bc in bcgraph.findall(".//boundary_condition"): bc_dict = bc.attrib if bc_dict['type'] == 'transmural pressures': @@ -257,10 +258,10 @@ def ReadFromXML(self, xmlBcpath, xsdBcpath=None): for entity,meshlist in self.NetworkMesh.Entities.iteritems(): if entity.Id == ent_dict['id']: for mesh in meshlist: - self.PressureValues[mesh.Id] = self.TimePressure + self.PressureValues[mesh.Id] = self.TimePressure if data.tag == "pressure": for pressure in data.findall(".//scalar"): - pressure_v = float(pressure.text) + pressure_v = float(pressure.text) for entities in bc.findall(".//entities"): if bc_dict['id'] == id: for ent in entities.findall(".//entity"): @@ -272,57 +273,61 @@ def ReadFromXML(self, xmlBcpath, xsdBcpath=None): raise EntityDuplicateError(entity) else: self.PressureValues[mesh.Id] = pressure_v - if bc_dict['type'] == 'input pressure': + if bc_dict['type'] == 'input pressure': id = bc_dict['id'] for param in bc.findall(".//parameters"): for data in param: - if data.tag == "pressure": + if data.tag == "pressure": for pressure in data.findall(".//scalar"): pressure_vp = float(pressure.text) self.InP = pressure_vp - for entities in bc.findall(".//entities"): + for entities in bc.findall(".//entities"): for ent in entities.findall(".//entity"): ent_dict = ent.attrib - ent_venp = ent_dict['id'] + ent_venp = ent_dict['id'] for entities in self.NetworkMesh.Entities.iterkeys(): - if ent_venp == entities.Id: - elNodesList = [] + if ent_venp == entities.Id: + elNodesList = [] for el in self.NetworkMesh.Entities[entities]: - elNodesList.append(el.NodeIds[1]) - self.NodeIn = min(elNodesList) + elNodesList.append(el.NodeIds[1]) + self.NodeIn = min(elNodesList) for el in self.NetworkMesh.Elements: - if el.NodeIds[1] == self.NodeIn: - self.elementIn = el - - if bc_dict['type'] == 'outflow pressure': + if el.NodeIds[1] == self.NodeIn: + self.elementIn = el + + if bc_dict['type'] == 'outflow pressure': id = bc_dict['id'] for param in bc.findall(".//parameters"): for data in param: - if data.tag == "pressure": + if data.tag == "pressure": for pressure in data.findall(".//scalar"): pressure_vp = float(pressure.text) - self.OutP = pressure_vp - for entities in bc.findall(".//entities"): + self.OutP = True + for entities in bc.findall(".//entities"): for ent in entities.findall(".//entity"): ent_dict = ent.attrib - ent_venp = ent_dict['id'] + ent_venp = ent_dict['id'] for entities in self.NetworkMesh.Entities.iterkeys(): - if ent_venp == entities.Id: - elNodesList = [] + if ent_venp == entities.Id: + elNodesList = [] for el in self.NetworkMesh.Entities[entities]: - elNodesList.append(el.NodeIds[1]) - self.NodeOut = max(elNodesList) + elNodesList.append(el.NodeIds[1]) + self.NodeOut= max(elNodesList) for el in self.NetworkMesh.Elements: - if el.NodeIds[1] == self.NodeOut: - self.elementOut = el - + if el.NodeIds[1] == self.NodeOut: + self.PressureOut[el] = {'node': self.NodeOut, 'value': pressure_vp} + + + + + if bc_dict['type'] == 'inflow': - - for entities in bc.findall(".//entities"): + + for entities in bc.findall(".//entities"): for ent in entities.findall(".//entity"): ent_dict = ent.attrib - ent_flow = ent_dict['id'] - node_flow = ent_dict['node_id'] + ent_flow = ent_dict['id'] + node_flow = ent_dict['node_id'] for entities in self.NetworkMesh.Entities.iterkeys(): if ent_flow == entities.Id: for el in self.NetworkMesh.Entities[entities]: @@ -335,24 +340,24 @@ def ReadFromXML(self, xmlBcpath, xsdBcpath=None): for data in param: if data.tag == "A0": for A0 in data.findall(".//scalar"): - + self.InFlows[el]['A0'] = float(A0.text) - + if data.tag == "fourier_coeffs": f_dict = data.attrib n = int(f_dict['n']) m = int(f_dict['m']) f_coeff = zeros((n,m)) - for fourier_coeffs in data.findall(".//matrix_nxm"): + for fourier_coeffs in data.findall(".//matrix_nxm"): f_coeff = array(fourier_coeffs.text.split(), dtype = float).reshape(m,n) self.InFlows[el]['f_coeff'] = f_coeff - + if data.tag == "signal": for values in data.findall(".//values"): self.InFlows[el]['signal'] = values.text.split() - - - + + + class Error(Exception): ''' A base class for exceptions defined in this module. @@ -367,27 +372,27 @@ def __init__(self,xmlschema): print "Error, Invalid Boundary Condition Xml File." print xmlschema.error_log sys.exit() - + class NoXSDWarning(Error): ''' Exception raised if no xsd file is provided. ''' def __init__(self): print "Warning, XML schema file was not provided.\nBoundary Conditions Xml file can not be validated." - + class WrongXSDPathError(Error): ''' Exception raised if a wrong xsd path is provided. ''' def __init__(self): print "Warning, Xml schema file not found.\nBoundary Conditions Xml file can not be validated." - + class LXMLError(Error): ''' Exception raised if lxml package is not installed. ''' def __init__(self): - print "Warning, Lxml package was not provided.\nBoundary Conditions Xml file can not be validated." + print "Warning, Lxml package was not provided.\nBoundary Conditions Xml file can not be validated." class XMLIdError(Error): ''' @@ -402,4 +407,4 @@ class EntityDuplicateError(Error): one boundary condition of the same type. ''' def __init__(self, entity): - print "This entity (", entity.Id, ") is already specified for another boundary condition of the same type.\nCheck your Boundary Conditions XML File" \ No newline at end of file + print "This entity (", entity.Id, ") is already specified for another boundary condition of the same type.\nCheck your Boundary Conditions XML File" diff --git a/Solver.py b/Solver.py index 5c9081a..ec15c79 100644 --- a/Solver.py +++ b/Solver.py @@ -9,8 +9,8 @@ ## Copyright (c) Simone Manini, Luca Antiga. All rights reserved. ## See LICENCE file for details. -## This software is distributed WITHOUT ANY WARRANTY; without even -## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR ## PURPOSE. See the above copyright notices for more information. ## Developed with support from the EC FP7/2007-2013: ARCH, Project n. 224390 @@ -30,7 +30,7 @@ class Solver(object): This is a general Solver Class. It doesn't provide any solver methods. This class provides only parameters setting methods ''' - + def __init__(self): ''' Constructor @@ -45,66 +45,66 @@ def __init__(self): self.Period = None self.CardiacFreq = None self.Cycles = None - self.SimulationDays = [] # Days' list for adaptation - self.NumberOfIncrements = None + self.SimulationDays = [] # Days' list for adaptation + self.NumberOfIncrements = None self.IncrementNumber = 1 # increment number self.EndIncrementTime = 0.0 # end increment time self.nltol = float(1e-3) # nonLinear convergence criterium - self.convergence = float(1e-4) # steady convergence limit for steady pre runs + self.convergence = float(1e-4) # steady convergence limit for steady pre runs self.Flow = None self.PrescribedPressures = None self.LinearZeroOrderGlobalMatrix = None self.LinearFirstOrderGlobalMatrix = None self.LinearSecondOrderGlobalMatrix = None - + def SetNetworkMesh(self,networkMesh): ''' Setting NetworkMesh ''' self.NetworkMesh = networkMesh - + def SetEvaluator(self, evaluator): ''' Setting Evaluator ''' self.Evaluator = evaluator - + def SetSimulationContext(self, simulationContext): ''' Setting SimulationContext ''' self.SimulationContext = simulationContext - + def SetBoundaryConditions(self, boundaryConditions): ''' Setting BoundaryConditions ''' self.BoundaryConditions = boundaryConditions - + def SetSteadyFlow(self): ''' Setting Steady Flow ''' self.steady = True - + def SetPulseFlow(self): ''' Setting Pulsatile Flow ''' self.steady = False - + def SetNonLinearTolerance(self, nltol): ''' Setting Non Linear Tolerance Value ''' self.nltol = float(nltol) - + def SetSteadyConvergenceLimit(self, convergence): ''' Setting convergence limit for steady pre-run simulations ''' self.convergence = convergence - + class SolverFirstTrapezoid(Solver): ''' This class provide a method to solving the system with "First Order Trapezium Method" @@ -115,13 +115,13 @@ def __init__(self): Constructor ''' Solver.__init__(self) - + def Solve(self): ''' This method builds System Matrix and gets Solution ''' if self.SimulationContext.Id != self.NetworkMesh.Id: - raise self.SimulationContext.XMLIdError() + raise self.SimulationContext.XMLIdError() try: self.TimeStep = self.SimulationContext.Context['timestep'] self.SquareTimeStep = self.TimeStep*self.TimeStep @@ -140,27 +140,27 @@ def Solve(self): except KeyError: print "Error, Please set cycles number in Simulation Context XML File" raise - + history = [] assembler = Assembler() assembler.SetNetworkMesh(self.NetworkMesh) assembler.SetBoundaryConditions(self.BoundaryConditions) info = {'dofmap':assembler.DofMap,'solution':None,'incrementNumber':self.IncrementNumber,'history':history} self.Evaluator.SetInfo(info) - + self.PrescribedPressures = assembler.AssembleBoundaryConditions(self.SimulationContext) - + self.LinearZeroOrderGlobalMatrix, self.LinearFirstOrderGlobalMatrix, self.LinearSecondOrderGlobalMatrix = \ assembler.AssembleInit(self.SimulationContext, self.Evaluator) - + self.ZeroOrderGlobalMatrix = assembler.ZeroOrderGlobalMatrix self.FirstOrderGlobalMatrix = assembler.FirstOrderGlobalMatrix - self.SecondOrderGlobalMatrix = assembler.SecondOrderGlobalMatrix + self.SecondOrderGlobalMatrix = assembler.SecondOrderGlobalMatrix - NumberOfGlobalDofs = assembler.GetNumberOfGlobalDofs() # number of dofs - self.UnknownPressures = arange(0,NumberOfGlobalDofs).reshape(NumberOfGlobalDofs,1) # unknown pressures + NumberOfGlobalDofs = assembler.GetNumberOfGlobalDofs() # number of dofs + self.UnknownPressures = arange(0,NumberOfGlobalDofs).reshape(NumberOfGlobalDofs,1) # unknown pressures self.UnknownPressures = delete(self.UnknownPressures, s_[self.PrescribedPressures[:,0]], axis=0) - PressuresMatrix = zeros((NumberOfGlobalDofs, self.NumberOfIncrements)) + PressuresMatrix = zeros((NumberOfGlobalDofs, self.NumberOfIncrements)) self.p = zeros((NumberOfGlobalDofs,1)) self.pt = zeros((NumberOfGlobalDofs,1)) self.ptt = zeros((NumberOfGlobalDofs,1)) @@ -182,16 +182,16 @@ def Solve(self): nonLinear = True break - while self.IncrementNumber<=self.NumberOfIncrements: + while self.IncrementNumber<=self.NumberOfIncrements: icc = (self.IncrementNumber%self.TimeStepFreq) if icc == 0: icc = self.TimeStepFreq - + #for flow in self.BoundaryConditions.elementFlow: for el in self.BoundaryConditions.elementFlow: if self.steady == True: self.Flow = assembler.BoundaryConditions.GetSteadyFlow(el, self.TimeStep,icc*self.TimeStep) - else: + else: self.Flow = assembler.BoundaryConditions.GetTimeFlow(el, icc*self.TimeStep) self.fe[assembler.FlowDof[el.Id]]= self.Flow @@ -202,29 +202,29 @@ def Solve(self): sumvbk[:,:] = self.sumv[:,:] counter = 0 while True: - #Build the algebric equation system for the increment + #Build the algebric equation system for the increment SystemMatrix = (2.0/self.TimeStep)*self.SecondOrderGlobalMatrix + self.FirstOrderGlobalMatrix + (self.TimeStep/2.0)*self.ZeroOrderGlobalMatrix #system matrix - RightVector = self.fe + (2.0/self.TimeStep)*dot(self.SecondOrderGlobalMatrix,(self.pt)) + dot(self.SecondOrderGlobalMatrix,(self.dpt)) - dot(self.ZeroOrderGlobalMatrix,(self.sumv))-(self.TimeStep/2.0)*dot(self.ZeroOrderGlobalMatrix,(self.pt)) # right hand side vector - #The reduced (partioned) system of equations is generated. + RightVector = self.fe + (2.0/self.TimeStep)*dot(self.SecondOrderGlobalMatrix,(self.pt)) + dot(self.SecondOrderGlobalMatrix,(self.dpt)) - dot(self.ZeroOrderGlobalMatrix,(self.sumv))-(self.TimeStep/2.0)*dot(self.ZeroOrderGlobalMatrix,(self.pt)) # right hand side vector + #The reduced (partioned) system of equations is generated. RightVector[:,:] = RightVector[:,:] - dot(SystemMatrix[:,self.PrescribedPressures[:,0]],self.PrescribedPressures[:,1:]) SystemMatrix = SystemMatrix[:,s_[self.UnknownPressures[:,0]]] - if SystemMatrix.shape[0]> 0.0: + if SystemMatrix.shape[0]> 0.0: SystemMatrix = SystemMatrix[s_[self.UnknownPressures[:,0]],:] RightVector = RightVector[s_[self.UnknownPressures[:,0]],:] #Unknown nodal point values are solved from this system. # Prescribed nodal values are inserted in the solution vector. - Solution = solve(SystemMatrix,RightVector) # solutions, unknown pressures - self.p[self.UnknownPressures,0] = Solution[:,:] + Solution = solve(SystemMatrix,RightVector) # solutions, unknown pressures + self.p[self.UnknownPressures,0] = Solution[:,:] self.p[self.PrescribedPressures[:,0],0] = self.PrescribedPressures[:,1] #Calculating derivatives. #Calculating internal nodal flow values. self.dp = dot((2.0/self.TimeStep),(self.p-self.pt)) - self.dpt self.ddp = dot((4.0/self.SquareTimeStep),(self.p-self.pt)) - dot((4.0/self.TimeStep),self.dpt) -self.ddpt self.sumv = sumvbk + dot((self.TimeStep/2.0),(self.pt+self.p)) - self.fi = dot(self.SecondOrderGlobalMatrix,(self.dp)) + dot(self.FirstOrderGlobalMatrix,(self.p)) + dot(self.ZeroOrderGlobalMatrix,(self.sumv)) + self.fi = dot(self.SecondOrderGlobalMatrix,(self.dp)) + dot(self.FirstOrderGlobalMatrix,(self.p)) + dot(self.ZeroOrderGlobalMatrix,(self.sumv)) if not nonLinear : break - + if self.pi == None: self.pi = zeros((NumberOfGlobalDofs,1)) self.pi[:,:] = self.pt[:,:] @@ -234,16 +234,16 @@ def Solve(self): if den < 1e-12: den = 1.0 nlerr = norm(self.p-self.pi,Inf) / den - + info = {'dofmap':assembler.DofMap,'solution':[self.p, self.pt, self.ptt],'incrementNumber':self.IncrementNumber,'history':history} self.Evaluator.SetInfo(info) assembler.Assemble(self.SimulationContext, self.Evaluator, self.LinearZeroOrderGlobalMatrix, self.LinearFirstOrderGlobalMatrix, self.LinearSecondOrderGlobalMatrix) self.ZeroOrderGlobalMatrix = assembler.ZeroOrderGlobalMatrix self.FirstOrderGlobalMatrix = assembler.FirstOrderGlobalMatrix - self.SecondOrderGlobalMatrix = assembler.SecondOrderGlobalMatrix - - #Dynamic nonlinear relaxing coefficient + self.SecondOrderGlobalMatrix = assembler.SecondOrderGlobalMatrix + + #Dynamic nonlinear relaxing coefficient if counter == 100: print "relaxing..." print nlerr, nltol, CoeffRelax @@ -254,25 +254,25 @@ def Solve(self): nltol *= 0.95 if nlerr < nltol: nltol = self.nltol - counter = 0 + counter = 0 break counter+=1 self.pi[:,:] = self.p[:,:] - + self.ptt[:,:] = self.pt[:,:] self.pt[:,:] = self.p[:,:] self.dpt[:,:] = self.dp[:,:] self.ddpt[:,:] = self.ddp[:,:] self.fet[:,:] = self.fe[:,:] - self.fit[:,:] = self.fi[:,:] - PressuresMatrix[:,(self.IncrementNumber-1)] = self.p[:,0] + self.fit[:,:] = self.fi[:,:] + PressuresMatrix[:,(self.IncrementNumber-1)] = self.p[:,0] history.insert(0,self.IncrementNumber) history = history[:3] - + if self.steady == True: self.MinimumIncrementNumber = 0.01* self.NumberOfIncrements if norm(self.fi-self.fe,Inf) self.MinimumIncrementNumber: - self.IncrementNumber = self.NumberOfIncrements + self.IncrementNumber = self.NumberOfIncrements else: pass @@ -281,17 +281,17 @@ def Solve(self): if self.IncrementNumber==ceil(0.25*self.NumberOfIncrements): print "->25%" if self.IncrementNumber==ceil(0.5*self.NumberOfIncrements): - print "->50%" + print "->50%" if self.IncrementNumber==ceil(0.70*self.NumberOfIncrements): - print "->70%" + print "->70%" if self.IncrementNumber==ceil(0.90*self.NumberOfIncrements): - print "->90%" + print "->90%" if self.IncrementNumber==ceil(0.99*self.NumberOfIncrements): - print "->99%" - + print "->99%" + self.IncrementNumber = self.IncrementNumber+1 - self.EndIncrementTime = self.EndIncrementTime + self.TimeStep # increment - info = {'dofmap':assembler.DofMap,'solution':[self.p, self.pt, self.ptt],'incrementNumber':self.IncrementNumber,'history':history,'allSolution':PressuresMatrix} + self.EndIncrementTime = self.EndIncrementTime + self.TimeStep # increment + info = {'dofmap':assembler.DofMap,'solution':[self.p, self.pt, self.ptt],'incrementNumber':self.IncrementNumber,'history':history,'allSolution':PressuresMatrix} self.Evaluator.SetInfo(info) self.Solutions = PressuresMatrix - return PressuresMatrix \ No newline at end of file + return PressuresMatrix