PageRenderTime 108ms CodeModel.GetById 20ms app.highlight 82ms RepoModel.GetById 1ms app.codeStats 0ms

/FiniteElasticity/SCRATCH/QuadraticEllipsoidCosta/src/QuadraticEllipsoidCostaExample.f90

http://github.com/xyan075/examples
FORTRAN Modern | 629 lines | 430 code | 76 blank | 123 comment | 0 complexity | 7bbbfcf1652c989237790595f2ec4938 MD5 | raw file
  1!> \file
  2!> $Id: QuadraticEllipsoidExample.f90 1714 2010-11-10 15:43:15Z jackisawesome $
  3!> \author Chris Bradley
  4!> \brief This is an example program to solve a finite elasticity equation using openCMISS calls.
  5!>
  6!> \section LICENSE
  7!>
  8!> Version: MPL 1.1/GPL 2.0/LGPL 2.1
  9!>
 10!> The contents of this file are subject to the Mozilla Public License
 11!> Version 1.1 (the "License"); you may not use this file except in
 12!> compliance with the License. You may obtain a copy of the License at
 13!> http://www.mozilla.org/MPL/
 14!>
 15!> Software distributed under the License is distributed on an "AS IS"
 16!> basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
 17!> License for the specific language governing rights and limitations
 18!> under the License.
 19!>
 20!> The Original Code is openCMISS
 21!>
 22!> The Initial Developer of the Original Code is University of Auckland,
 23!> Auckland, New Zealand and University of Oxford, Oxford, United
 24!> Kingdom. Portions created by the University of Auckland and University
 25!> of Oxford are Copyright (C) 2007 by the University of Auckland and
 26!> the University of Oxford. All Rights Reserved.
 27!>
 28!> Contributor(s): Jack Lee
 29!>
 30!> Alternatively, the contents of this file may be used under the terms of
 31!> either the GNU General Public License Version 2 or later (the "GPL"), or
 32!> the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
 33!> in which case the provisions of the GPL or the LGPL are applicable instead
 34!> of those above. If you wish to allow use of your version of this file only
 35!> under the terms of either the GPL or the LGPL, and not to allow others to
 36!> use your version of this file under the terms of the MPL, indicate your
 37!> decision by deleting the provisions above and replace them with the notice
 38!> and other provisions required by the GPL or the LGPL. If you do not delete
 39!> the provisions above, a recipient may use your version of this file under
 40!> the terms of any one of the MPL, the GPL or the LGPL.
 41!>
 42
 43!> \example FiniteElasticity/UniAxialExtension/src/UniAxialExtensionExample.f90
 44!! Example program to solve a finite elasticity equation using openCMISS calls.
 45!! \par Latest Builds:
 46!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/FiniteElasticity/UniAxialExtension/build-intel'>Linux Intel Build</a>
 47!! \li <a href='http://autotest.bioeng.auckland.ac.nz/opencmiss-build/logs_x86_64-linux/FiniteElasticity/UniAxialExtension/build-gnu'>Linux GNU Build</a>
 48!<
 49
 50!> Main program
 51PROGRAM QUADRATICELLIPSOIDCOSTAEXAMPLE
 52
 53  USE OPENCMISS
 54  USE MPI
 55
 56#ifdef WIN32
 57  USE IFQWIN
 58#endif
 59
 60  IMPLICIT NONE
 61
 62  !Test program parameters
 63
 64  REAL(CMISSDP), PARAMETER :: PI=3.14159_CMISSDP
 65  REAL(CMISSDP), PARAMETER :: LONG_AXIS=2.0_CMISSDP
 66  REAL(CMISSDP), PARAMETER :: SHORT_AXIS=1.0_CMISSDP
 67  REAL(CMISSDP), PARAMETER :: WALL_THICKNESS=0.5_CMISSDP
 68  REAL(CMISSDP), PARAMETER :: CUTOFF_ANGLE=1.5708_CMISSDP
 69  REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_ENDO=1.73205_CMISSDP !Slope of fibres in base endocardium = 60 degrees
 70  REAL(CMISSDP), PARAMETER :: FIBRE_SLOPE_EPI=-1.73205_CMISSDP !Slope change of fibres from 60 to -60 degrees in transmural direction
 71  REAL(CMISSDP), PARAMETER :: SHEET_SLOPE_BASE_ENDO=1.0_CMISSDP !Slope of sheet at base endocardium
 72  REAL(CMISSDP), DIMENSION(1:7) :: COSTA_PARAMS =  (/ 0.2, 30.0, 12.0, 14.0, 14.0, 10.0, 18.0 /) ! a bff bfs bfn bss bsn bnn
 73  REAL(CMISSDP), PARAMETER :: INNER_PRESSURE=1.0_CMISSDP  !Positive is compressive
 74  REAL(CMISSDP), PARAMETER :: OUTER_PRESSURE=0.0_CMISSDP  !Positive is compressive
 75
 76  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalXElements=4  ! X ==NUMBER_GLOBAL_CIRCUMFERENTIAL_ELEMENTS
 77  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalYElements=4  ! Y ==NUMBER_GLOBAL_LONGITUDINAL_ELEMENTS
 78  INTEGER(CMISSIntg), PARAMETER :: NumberGlobalZElements=1  ! Z ==NUMBER_GLOBAL_TRANSMURAL_ELEMENTS
 79  INTEGER(CMISSIntg) :: NumberOfDomains
 80
 81  INTEGER(CMISSIntg), PARAMETER :: CoordinateSystemUserNumber=1
 82  INTEGER(CMISSIntg), PARAMETER :: NumberOfSpatialCoordinates=3
 83  INTEGER(CMISSIntg), PARAMETER :: RegionUserNumber=1
 84  INTEGER(CMISSIntg), PARAMETER :: QuadraticBasisUserNumber=1
 85  INTEGER(CMISSIntg), PARAMETER :: QuadraticCollapsedBasisUserNumber=2
 86  INTEGER(CMISSIntg), PARAMETER :: LinearBasisUserNumber=3
 87  INTEGER(CMISSIntg), PARAMETER :: LinearCollapsedBasisUserNumber=4
 88  INTEGER(CMISSIntg), PARAMETER :: MeshUserNumber=1
 89  INTEGER(CMISSIntg), PARAMETER :: GeneratedMeshUserNumber=2
 90  INTEGER(CMISSIntg), PARAMETER :: DecompositionUserNumber=1
 91  INTEGER(CMISSIntg), PARAMETER :: DerivativeUserNumber=1
 92
 93  INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshDimensions=3
 94  INTEGER(CMISSIntg), PARAMETER :: NumberOfXiCoordinates=3
 95  INTEGER(CMISSIntg), PARAMETER :: NumberOfMeshComponents=2
 96  INTEGER(CMISSIntg), PARAMETER :: QuadraticMeshComponentNumber=1
 97  INTEGER(CMISSIntg), PARAMETER :: LinearMeshComponentNumber=2
 98  INTEGER(CMISSIntg), PARAMETER :: TotalNumberOfElements=1
 99
100  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryUserNumber=1
101  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfVariables=1
102  INTEGER(CMISSIntg), PARAMETER :: FieldGeometryNumberOfComponents=3
103
104  INTEGER(CMISSIntg), PARAMETER :: FieldFibreUserNumber=2
105  INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfVariables=1
106  INTEGER(CMISSIntg), PARAMETER :: FieldFibreNumberOfComponents=3
107
108  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialUserNumber=3
109  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfVariables=1
110  INTEGER(CMISSIntg), PARAMETER :: FieldMaterialNumberOfComponents=7
111
112  INTEGER(CMISSIntg), PARAMETER :: FieldDependentUserNumber=4
113  INTEGER(CMISSIntg), PARAMETER :: FieldDependentNumberOfVariables=2
114  INTEGER(CMISSIntg), PARAMETER :: FieldDependentNumberOfComponents=4
115
116  INTEGER(CMISSIntg), PARAMETER :: EquationSetUserNumber=1
117  INTEGER(CMISSIntg), PARAMETER :: EquationsSetFieldUserNumber=13
118  INTEGER(CMISSIntg), PARAMETER :: ProblemUserNumber=1
119
120  !Program types
121
122
123  !Program variables
124
125  INTEGER(CMISSIntg) :: MPI_IERROR
126  INTEGER(CMISSIntg) :: EquationsSetIndex
127  INTEGER(CMISSIntg) :: NumberOfComputationalNodes,ComputationalNodeNumber
128  REAL(CMISSDP) :: FibreFieldAngle(3)
129  REAL(CMISSDP) :: nu,theta,omega,XI3,XI3delta,XI2delta, zero
130  INTEGER(CMISSIntg) ::i,j,k,time_step,component_idx,node_idx,TOTAL_NUMBER_NODES_XI(3)
131  CHARACTER(LEN=50) :: fileName,fileNumber
132  !For grabbing surfaces
133  INTEGER(CMISSIntg) :: InnerNormalXi,OuterNormalXi,TopNormalXi
134  INTEGER(CMISSIntg), ALLOCATABLE :: InnerSurfaceNodes(:)
135  INTEGER(CMISSIntg), ALLOCATABLE :: OuterSurfaceNodes(:)
136  INTEGER(CMISSIntg), ALLOCATABLE :: TopSurfaceNodes(:)
137  INTEGER(CMISSIntg), ALLOCATABLE :: G(:)
138  INTEGER(CMISSIntg) :: NotCornerNode, CornerNode, NumberOfCornerNodes, gn, CorrectNodeNumber,TotalNumberOfNodes
139  INTEGER(CMISSIntg) :: NN,NODE,NodeDomain
140  REAL(CMISSDP) :: XCoord,YCoord,ZCoord
141  LOGICAL :: X_FIXED,Y_FIXED,X_OKAY,Y_OKAY
142
143  !CMISS variables
144
145  TYPE(CMISSBasisType) :: QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis
146  TYPE(CMISSBoundaryConditionsType) :: BoundaryConditions
147  TYPE(CMISSCoordinateSystemType) :: CoordinateSystem, WorldCoordinateSystem
148  TYPE(CMISSMeshType) :: Mesh
149  TYPE(CMISSGeneratedMeshType) :: GeneratedMesh
150  TYPE(CMISSDecompositionType) :: Decomposition
151  TYPE(CMISSEquationsType) :: Equations
152  TYPE(CMISSEquationsSetType) :: EquationsSet
153  TYPE(CMISSFieldType) :: GeometricField,FibreField,MaterialField,DependentField,EquationsSetField
154  TYPE(CMISSFieldsType) :: Fields
155  TYPE(CMISSProblemType) :: Problem
156  TYPE(CMISSRegionType) :: Region,WorldRegion
157  TYPE(CMISSSolverType) :: Solver,LinearSolver
158  TYPE(CMISSSolverEquationsType) :: SolverEquations
159
160#ifdef WIN32
161  !Quickwin type
162  LOGICAL :: QUICKWIN_STATUS=.FALSE.
163  TYPE(WINDOWCONFIG) :: QUICKWIN_WINDOW_CONFIG
164#endif
165
166  !Generic CMISS variables
167  INTEGER(CMISSIntg) :: Err
168
169#ifdef WIN32
170  !Initialise QuickWin
171  QUICKWIN_WINDOW_CONFIG%TITLE="General Output" !Window title
172  QUICKWIN_WINDOW_CONFIG%NUMTEXTROWS=-1 !Max possible number of rows
173  QUICKWIN_WINDOW_CONFIG%MODE=QWIN$SCROLLDOWN
174  !Set the window parameters
175  QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
176  !If attempt fails set with system estimated values
177  IF(.NOT.QUICKWIN_STATUS) QUICKWIN_STATUS=SETWINDOWCONFIG(QUICKWIN_WINDOW_CONFIG)
178#endif
179
180  !Intialise cmiss
181  CALL CMISSInitialise(WorldCoordinateSystem,WorldRegion,Err)
182
183  CALL CMISSErrorHandlingModeSet(CMISS_ERRORS_TRAP_ERROR,Err)
184
185  WRITE(*,'(A)') "Program starting."
186
187  !Set all diganostic levels on for testing
188  CALL CMISSDiagnosticsSetOn(CMISS_FROM_DIAG_TYPE,(/1,2,3,4,5/),"Diagnostics",(/"PROBLEM_FINITE_ELEMENT_CALCULATE"/),Err)
189
190  !Get the number of computational nodes and this computational node number
191  CALL CMISSComputationalNumberOfNodesGet(NumberOfComputationalNodes,Err)
192  CALL CMISSComputationalNodeNumberGet(ComputationalNodeNumber,Err)
193
194  !Broadcast the number of elements in the X,Y and Z directions and the number of partitions to the other computational nodes
195!   CALL MPI_BCAST(NumberGlobalXElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
196!   CALL MPI_BCAST(NumberGlobalYElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
197!   CALL MPI_BCAST(NumberGlobalZElements,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
198!   CALL MPI_BCAST(NumberOfDomains,1,MPI_INTEGER,0,MPI_COMM_WORLD,MPI_IERROR)
199  NumberOfDomains=NumberOfComputationalNodes
200
201  !Create a CS - default is 3D rectangular cartesian CS with 0,0,0 as origin
202  CALL CMISSCoordinateSystem_Initialise(CoordinateSystem,Err)
203  CALL CMISSCoordinateSystem_CreateStart(CoordinateSystemUserNumber,CoordinateSystem,Err)
204  CALL CMISSCoordinateSystem_TypeSet(CoordinateSystem,CMISS_COORDINATE_RECTANGULAR_CARTESIAN_TYPE,Err)
205  CALL CMISSCoordinateSystem_DimensionSet(CoordinateSystem,NumberOfSpatialCoordinates,Err)
206  CALL CMISSCoordinateSystem_OriginSet(CoordinateSystem,(/0.0_CMISSDP,0.0_CMISSDP,0.0_CMISSDP/),Err)
207  CALL CMISSCoordinateSystem_CreateFinish(CoordinateSystem,Err)
208
209  !Create a region and assign the CS to the region
210  CALL CMISSRegion_Initialise(Region,Err)
211  CALL CMISSRegion_CreateStart(RegionUserNumber,WorldRegion,Region,Err)
212  CALL CMISSRegion_CoordinateSystemSet(Region,CoordinateSystem,Err)
213  CALL CMISSRegion_CreateFinish(Region,Err)
214
215  !Define basis functions - tri-linear Lagrange and tri-Quadratic Lagrange, each with collapsed variant
216    !Quadratic Basis
217  CALL CMISSBasis_Initialise(QuadraticBasis,Err)
218  CALL CMISSBasis_CreateStart(QuadraticBasisUserNumber,QuadraticBasis,Err)
219  CALL CMISSBasis_InterpolationXiSet(QuadraticBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
220    & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
221  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticBasis, &
222    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
223  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticBasis,.true.,Err) !Have to do this
224  CALL CMISSBasis_CreateFinish(QuadraticBasis,Err)
225
226    !Collapsed Quadratic Basis
227  CALL CMISSBasis_Initialise(QuadraticCollapsedBasis,Err)
228  CALL CMISSBasis_CreateStart(QuadraticCollapsedBasisUserNumber,QuadraticCollapsedBasis,Err)
229  CALL CMISSBasis_TypeSet(QuadraticCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
230  CALL CMISSBasis_NumberOfXiSet(QuadraticCollapsedBasis,NumberOfXiCoordinates,Err)
231  CALL CMISSBasis_InterpolationXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION, &
232       & CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION,CMISS_BASIS_QUADRATIC_LAGRANGE_INTERPOLATION/),Err)
233  CALL CMISSBasis_CollapsedXiSet(QuadraticCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED, &
234       & CMISS_BASIS_COLLAPSED_AT_XI0,CMISS_BASIS_NOT_COLLAPSED/),Err)
235  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(QuadraticCollapsedBasis, &
236       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
237  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(QuadraticCollapsedBasis,.true.,Err) !Have to do this
238  CALL CMISSBasis_CreateFinish(QuadraticCollapsedBasis,Err)
239
240    !Linear Basis
241  CALL CMISSBasis_Initialise(LinearBasis,Err)
242  CALL CMISSBasis_CreateStart(LinearBasisUserNumber,LinearBasis,Err)
243  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearBasis, &
244    & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
245  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
246  CALL CMISSBasis_CreateFinish(LinearBasis,Err)
247
248    !Collapsed Linear Basis
249  CALL CMISSBasis_Initialise(LinearCollapsedBasis,Err)
250  CALL CMISSBasis_CreateStart(LinearCollapsedBasisUserNumber,LinearCollapsedBasis,Err)
251  CALL CMISSBasis_TypeSet(LinearCollapsedBasis,CMISS_BASIS_LAGRANGE_HERMITE_TP_TYPE,Err)
252  CALL CMISSBasis_NumberOfXiSet(LinearCollapsedBasis,NumberOfXiCoordinates,Err)
253  CALL CMISSBasis_InterpolationXiSet(LinearCollapsedBasis,(/CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION, &
254       & CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION,CMISS_BASIS_LINEAR_LAGRANGE_INTERPOLATION/),Err)
255  CALL CMISSBasis_CollapsedXiSet(LinearCollapsedBasis,(/CMISS_BASIS_XI_COLLAPSED,CMISS_BASIS_COLLAPSED_AT_XI0, &
256    & CMISS_BASIS_NOT_COLLAPSED/),Err)
257  CALL CMISSBasis_QuadratureNumberOfGaussXiSet(LinearCollapsedBasis, &
258       & (/CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME,CMISS_BASIS_MID_QUADRATURE_SCHEME/),Err)
259  CALL CMISSBasis_QuadratureLocalFaceGaussEvaluateSet(LinearCollapsedBasis,.true.,Err) !Have to do this (unused) due to field_interp setup
260  CALL CMISSBasis_CreateFinish(LinearCollapsedBasis,Err)
261
262  !Start the creation of a generated ellipsoid mesh
263  CALL CMISSGeneratedMesh_Initialise(GeneratedMesh,Err)
264  CALL CMISSGeneratedMesh_CreateStart(GeneratedMeshUserNumber,Region,GeneratedMesh,Err)
265  !Set up an ellipsoid mesh
266  CALL CMISSGeneratedMesh_TypeSet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_MESH_TYPE,Err)
267   !Set the quadratic and linear bases
268  CALL CMISSGeneratedMesh_BasisSet(GeneratedMesh,[QuadraticBasis,QuadraticCollapsedBasis,LinearBasis,LinearCollapsedBasis],Err)
269  !Define the mesh on the region
270  CALL CMISSGeneratedMesh_ExtentSet(GeneratedMesh,(/LONG_AXIS,SHORT_AXIS,WALL_THICKNESS,CUTOFF_ANGLE/),Err)
271  CALL CMISSGeneratedMesh_NumberOfElementsSet(GeneratedMesh,(/NumberGlobalXElements,NumberGlobalYElements, &
272    & NumberGlobalZElements/),Err)
273
274  !Finish the creation of a generated mesh in the region
275  CALL CMISSMesh_Initialise(Mesh,Err)
276  CALL CMISSGeneratedMesh_CreateFinish(GeneratedMesh,MeshUserNumber,Mesh,Err)
277
278  !Create a decomposition
279  CALL CMISSDecomposition_Initialise(Decomposition,Err)
280  CALL CMISSDecomposition_CreateStart(DecompositionUserNumber,Mesh,Decomposition,Err)
281  CALL CMISSDecomposition_TypeSet(Decomposition,CMISS_DECOMPOSITION_CALCULATED_TYPE,Err)
282  CALL CMISSDecomposition_NumberOfDomainsSet(Decomposition,NumberOfDomains,Err)
283  CALL CMISSDecomposition_CalculateFacesSet(Decomposition,.TRUE.,Err)
284  CALL CMISSDecomposition_CreateFinish(Decomposition,Err)
285
286  !Create a field to put the geometry (default is geometry)
287  CALL CMISSField_Initialise(GeometricField,Err)
288  CALL CMISSField_CreateStart(FieldGeometryUserNumber,Region,GeometricField,Err)
289  CALL CMISSField_MeshDecompositionSet(GeometricField,Decomposition,Err)
290  CALL CMISSField_TypeSet(GeometricField,CMISS_FIELD_GEOMETRIC_TYPE,Err)
291  CALL CMISSField_NumberOfVariablesSet(GeometricField,FieldGeometryNumberOfVariables,Err)
292  CALL CMISSField_NumberOfComponentsSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,FieldGeometryNumberOfComponents,Err)
293  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
294  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
295  CALL CMISSField_ComponentMeshComponentSet(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
296  CALL CMISSField_CreateFinish(GeometricField,Err)
297
298  !Update the geometric field parameters
299  CALL CMISSGeneratedMesh_GeometricParametersCalculate(GeneratedMesh,GeometricField,Err)
300
301  !Create a fibre field and attach it to the geometric field
302  CALL CMISSField_Initialise(FibreField,Err)
303  CALL CMISSField_CreateStart(FieldFibreUserNumber,Region,FibreField,Err)
304  CALL CMISSField_TypeSet(FibreField,CMISS_FIELD_FIBRE_TYPE,Err)
305  CALL CMISSField_MeshDecompositionSet(FibreField,Decomposition,Err)
306  CALL CMISSField_GeometricFieldSet(FibreField,GeometricField,Err)
307  CALL CMISSField_NumberOfVariablesSet(FibreField,FieldFibreNumberOfVariables,Err)
308  CALL CMISSField_NumberOfComponentsSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,FieldFibreNumberOfComponents,Err)
309  CALL CMISSField_ComponentMeshComponentSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
310  CALL CMISSField_ComponentMeshComponentSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
311  CALL CMISSField_ComponentMeshComponentSet(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
312  CALL CMISSField_CreateFinish(FibreField,Err)
313
314  !Set Fibre directions (this block is parallel-untested)
315  node_idx=0
316  !This is valid only for quadratic basis functions
317  TOTAL_NUMBER_NODES_XI(1)=NumberGlobalXElements*2
318  TOTAL_NUMBER_NODES_XI(2)=NumberGlobalYElements*2+1
319  TOTAL_NUMBER_NODES_XI(3)=NumberGlobalZElements*2+1
320
321  !Map the correct node number (cn) to geometric node number (gn) G(gn)=cn
322  NumberOfCornerNodes=(NumberGlobalZElements+1)*(NumberGlobalYElements*NumberGlobalXElements+1)
323  TotalNumberOfNodes=(TOTAL_NUMBER_NODES_XI(3))*((TOTAL_NUMBER_NODES_XI(2)-1)*TOTAL_NUMBER_NODES_XI(1)+1)
324
325  ALLOCATE(G(TotalNumberOfNodes),STAT=Err)
326
327  CornerNode=0
328  !Numbering of not corner nodes starts where Corner nodes end
329  NotCornerNode=NumberOfCornerNodes
330  gn=0
331
332  DO k=1,TOTAL_NUMBER_NODES_XI(3)
333     gn=gn+1
334     IF (mod(k,2)==0) THEN
335        !Not corner node
336        NotCornerNode=NotCornerNode+1
337        G(gn)=NotCornerNode
338     ELSE
339        !Corner node
340        CornerNode=CornerNode+1
341        G(gn)=CornerNode
342     ENDIF
343     DO j=2,TOTAL_NUMBER_NODES_XI(2)
344        DO i=1,TOTAL_NUMBER_NODES_XI(1)
345            gn=gn+1
346            IF ((mod(i,2)==0).OR.(mod(j,2)==0).OR.(mod(k,2)==0)) THEN
347               !Not corner node
348               NotCornerNode=NotCornerNode+1
349               G(gn)=NotCornerNode
350            ELSE
351               !Corner node
352               CornerNode=CornerNode+1
353               G(gn)=CornerNode
354            ENDIF
355         ENDDO
356      ENDDO
357   ENDDO
358
359
360  XI2delta=(PI-CUTOFF_ANGLE)/(TOTAL_NUMBER_NODES_XI(2)-1)
361  XI3=0
362  XI3delta=(1.0)/(TOTAL_NUMBER_NODES_XI(3)-1)
363  zero=0
364  DO k=1, TOTAL_NUMBER_NODES_XI(3)
365    !Apex nodes
366    j=1
367    i=1
368    node_idx=node_idx+1
369    CorrectNodeNumber=G(node_idx)
370    CALL CMISSDecomposition_NodeDomainGet(Decomposition,CorrectNodeNumber,1,NodeDomain,Err)
371    IF(NodeDomain==ComputationalNodeNumber) THEN
372      FibreFieldAngle=(/zero,zero,zero/)
373      DO component_idx=1,FieldFibreNumberOfComponents
374        CALL CMISSField_ParameterSetUpdateNode(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
375          & DerivativeUserNumber, &
376          & CorrectNodeNumber,component_idx,FibreFieldAngle(component_idx),Err)
377      ENDDO
378    ENDIF
379    theta=atan((FIBRE_SLOPE_EPI-FIBRE_SLOPE_ENDO)*XI3+FIBRE_SLOPE_ENDO)
380    DO j=2, TOTAL_NUMBER_NODES_XI(2)
381      nu=PI-XI2delta*(j-1)
382      omega=PI/2+cos(2*nu)*atan(SHEET_SLOPE_BASE_ENDO)*(-2*XI3+1)
383      DO i=1, TOTAL_NUMBER_NODES_XI(1)
384        node_idx=node_idx+1
385        CorrectNodeNumber=G(node_idx)
386        CALL CMISSDecomposition_NodeDomainGet(Decomposition,CorrectNodeNumber,1,NodeDomain,Err)
387        IF(NodeDomain==ComputationalNodeNumber) THEN
388          FibreFieldAngle=(/theta,zero,omega/)
389          DO component_idx=1,FieldFibreNumberOfComponents
390            CALL CMISSField_ParameterSetUpdateNode(FibreField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1, &
391              & DerivativeUserNumber,CorrectNodeNumber,component_idx,FibreFieldAngle(component_idx),Err)
392          ENDDO
393        ENDIF
394      ENDDO
395    ENDDO
396    XI3=XI3+XI3delta
397  ENDDO
398
399  !Create a material field and attach it to the geometric field
400  CALL CMISSField_Initialise(MaterialField,Err)
401  CALL CMISSField_CreateStart(FieldMaterialUserNumber,Region,MaterialField,Err)
402  CALL CMISSField_TypeSet(MaterialField,CMISS_FIELD_MATERIAL_TYPE,Err)
403  CALL CMISSField_MeshDecompositionSet(MaterialField,Decomposition,Err)
404  CALL CMISSField_GeometricFieldSet(MaterialField,GeometricField,Err)
405  CALL CMISSField_NumberOfVariablesSet(MaterialField,FieldMaterialNumberOfVariables,Err)
406  CALL CMISSField_NumberOfComponentsSet(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,FieldMaterialNumberOfComponents,Err)
407!   CALL CMISSField_ComponentInterpolationSet(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,1,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
408!   CALL CMISSField_ComponentInterpolationSet(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,2,CMISS_FIELD_CONSTANT_INTERPOLATION,Err)
409  CALL CMISSField_CreateFinish(MaterialField,Err)
410
411  !Set Mooney-Rivlin constants c10 and c01 to 2.0 and 6.0 respectively.
412  !CALL CMISSField_ComponentValuesInitialise(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,2.0_CMISSDP,Err)
413  !CALL CMISSField_ComponentValuesInitialise(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,6.0_CMISSDP,Err)
414  !Set Costa material parameters
415  DO I=1,7
416   CALL CMISSField_ComponentValuesInitialise(MaterialField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,I, &
417     & COSTA_PARAMS(I),Err)
418  ENDDO
419  !Create the equations_set
420  CALL CMISSField_Initialise(EquationsSetField,Err)
421  CALL CMISSEquationsSet_Initialise(EquationsSet,Err)
422  CALL CMISSEquationsSet_CreateStart(EquationSetUserNumber,Region,FibreField,CMISS_EQUATIONS_SET_ELASTICITY_CLASS, &
423    & CMISS_EQUATIONS_SET_FINITE_ELASTICITY_TYPE,CMISS_EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE, &
424      & EquationsSetFieldUserNumber,&
425    & EquationsSetField,EquationsSet,Err)
426  ! CMISS_EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE
427  CALL CMISSEquationsSet_CreateFinish(EquationsSet,Err)
428
429  !Create the dependent field with 2 variables and 4 components (3 displacement, 1 pressure)
430  CALL CMISSField_Initialise(DependentField,Err)
431  CALL CMISSField_CreateStart(FieldDependentUserNumber,Region,DependentField,Err)
432  CALL CMISSField_TypeSet(DependentField,CMISS_FIELD_GENERAL_TYPE,Err)
433  CALL CMISSField_MeshDecompositionSet(DependentField,Decomposition,Err)
434  CALL CMISSField_GeometricFieldSet(DependentField,GeometricField,Err)
435  CALL CMISSField_DependentTypeSet(DependentField,CMISS_FIELD_DEPENDENT_TYPE,Err)
436  CALL CMISSField_NumberOfVariablesSet(DependentField,FieldDependentNumberOfVariables,Err)
437  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,FieldDependentNumberOfComponents,Err)
438  CALL CMISSField_NumberOfComponentsSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,FieldDependentNumberOfComponents,Err)
439  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
440  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
441  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
442  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err)
443  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,QuadraticMeshComponentNumber,Err)
444  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,2,QuadraticMeshComponentNumber,Err)
445  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,3,QuadraticMeshComponentNumber,Err)
446  CALL CMISSField_ComponentMeshComponentSet(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,4,LinearMeshComponentNumber,Err)
447  CALL CMISSField_ScalingTypeSet(DependentField,CMISS_FIELD_UNIT_SCALING,Err)
448  CALL CMISSField_CreateFinish(DependentField,Err)
449
450  CALL CMISSEquationsSet_DependentCreateStart(EquationsSet,FieldDependentUserNumber,DependentField,Err)
451  CALL CMISSEquationsSet_DependentCreateFinish(EquationsSet,Err)
452
453  CALL CMISSEquationsSet_MaterialsCreateStart(EquationsSet,FieldMaterialUserNumber,MaterialField,Err)
454  CALL CMISSEquationsSet_MaterialsCreateFinish(EquationsSet,Err)
455
456  !Create the equations set equations
457  CALL CMISSEquations_Initialise(Equations,Err)
458  CALL CMISSEquationsSet_EquationsCreateStart(EquationsSet,Equations,Err)
459  CALL CMISSEquations_SparsityTypeSet(Equations,CMISS_EQUATIONS_SPARSE_MATRICES,Err)
460  CALL CMISSEquations_OutputTypeSet(Equations,CMISS_EQUATIONS_NO_OUTPUT,Err)
461  CALL CMISSEquationsSet_EquationsCreateFinish(EquationsSet,Err)
462
463  !Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure
464  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
465    & 1,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,Err)
466  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
467    & 2,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,2,Err)
468  CALL CMISSField_ParametersToFieldParametersComponentCopy(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE, &
469    & 3,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,3,Err)
470  CALL CMISSField_ComponentValuesInitialise(DependentField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,4, &
471    & -14.0_CMISSDP, &
472    & Err)
473
474  !Define the problem
475  CALL CMISSProblem_Initialise(Problem,Err)
476  CALL CMISSProblem_CreateStart(ProblemUserNumber,Problem,Err)
477  CALL CMISSProblem_SpecificationSet(Problem,CMISS_PROBLEM_ELASTICITY_CLASS,CMISS_PROBLEM_FINITE_ELASTICITY_TYPE, &
478    & CMISS_PROBLEM_NO_SUBTYPE,Err)
479  CALL CMISSProblem_CreateFinish(Problem,Err)
480
481  !Create the problem control loop
482  CALL CMISSProblem_ControlLoopCreateStart(Problem,Err)
483  CALL CMISSProblem_ControlLoopCreateFinish(Problem,Err)
484
485  !Create the problem solvers
486  CALL CMISSSolver_Initialise(Solver,Err)
487  CALL CMISSSolver_Initialise(LinearSolver,Err)
488  CALL CMISSProblem_SolversCreateStart(Problem,Err)
489  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,1,Solver,Err)
490  CALL CMISSSolver_OutputTypeSet(Solver,CMISS_SOLVER_PROGRESS_OUTPUT,Err)
491  !CALL CMISSSolver_NewtonJacobianCalculationTypeSet(Solver,CMISS_SOLVER_NEWTON_JACOBIAN_FD_CALCULATED,Err)
492  CALL CMISSSolver_NewtonJacobianCalculationTypeSet(Solver,CMISS_SOLVER_NEWTON_JACOBIAN_EQUATIONS_CALCULATED,Err)
493  CALL CMISSSolver_NewtonLinearSolverGet(Solver,LinearSolver,Err)
494  CALL CMISSSolver_LinearTypeSet(LinearSolver,CMISS_SOLVER_LINEAR_DIRECT_SOLVE_TYPE,Err)
495  CALL CMISSProblem_SolversCreateFinish(Problem,Err)
496
497  !Create the problem solver equations
498  CALL CMISSSolver_Initialise(Solver,Err)
499  CALL CMISSSolverEquations_Initialise(SolverEquations,Err)
500  CALL CMISSProblem_SolverEquationsCreateStart(Problem,Err)
501  CALL CMISSProblem_SolverGet(Problem,CMISS_CONTROL_LOOP_NODE,1,Solver,Err)
502  CALL CMISSSolver_SolverEquationsGet(Solver,SolverEquations,Err)
503  CALL CMISSSolverEquations_SparsityTypeSet(SolverEquations,CMISS_SOLVER_SPARSE_MATRICES,Err)
504  CALL CMISSSolverEquations_EquationsSetAdd(SolverEquations,EquationsSet,EquationsSetIndex,Err)
505  CALL CMISSProblem_SolverEquationsCreateFinish(Problem,Err)
506
507  !Prescribe boundary conditions (absolute nodal parameters)
508  CALL CMISSBoundaryConditions_Initialise(BoundaryConditions,Err)
509  CALL CMISSSolverEquations_BoundaryConditionsCreateStart(SolverEquations,BoundaryConditions,Err)
510
511  !Grab the list of nodes on inner, outer and top surfaces
512  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_TOP_SURFACE,TopSurfaceNodes,TopNormalXi,Err)
513  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_INNER_SURFACE,InnerSurfaceNodes,InnerNormalXi,Err)
514  CALL CMISSGeneratedMesh_SurfaceGet(GeneratedMesh,CMISS_GENERATED_MESH_ELLIPSOID_OUTER_SURFACE,OuterSurfaceNodes,OuterNormalXi,Err)
515
516  ! ASSIGN BOUNDARY CONDITIONS ===============
517  !Fix base of the ellipsoid in z direction
518  DO NN=1,SIZE(TopSurfaceNodes,1)
519    NODE=TopSurfaceNodes(NN)
520    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
521    IF(NodeDomain==ComputationalNodeNumber) THEN
522      CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,3, &
523        & ZCoord, &
524        & Err)
525      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,3, &
526        & CMISS_BOUNDARY_CONDITION_FIXED,ZCoord,Err)
527    ENDIF
528  ENDDO
529
530  !Apply inner surface pressure
531  !NOTE: Surface pressure goes into pressure_values_set_type of the DELUDELN type
532  DO NN=1,SIZE(InnerSurfaceNodes,1)
533    NODE=InnerSurfaceNodes(NN)
534    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
535    IF(NodeDomain==ComputationalNodeNumber) THEN
536      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, &
537        & ABS(InnerNormalXi), &
538        & CMISS_BOUNDARY_CONDITION_PRESSURE,INNER_PRESSURE,Err)
539    ENDIF
540  ENDDO
541
542  !Apply outer surface pressure
543  DO NN=1,SIZE(OuterSurfaceNodes,1)
544    NODE=OuterSurfaceNodes(NN)
545    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
546    IF(NodeDomain==ComputationalNodeNumber) THEN
547      CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE,1,1,NODE, &
548        & ABS(OuterNormalXi), &
549        & CMISS_BOUNDARY_CONDITION_PRESSURE,OUTER_PRESSURE,Err)
550    ENDIF
551  ENDDO
552
553  !Fix more nodes at the base to stop free body motion
554  X_FIXED=.FALSE.
555  Y_FIXED=.FALSE.
556  DO NN=1,SIZE(TopSurfaceNodes,1)
557    NODE=TopSurfaceNodes(NN)
558    CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
559    IF(NodeDomain==ComputationalNodeNumber) THEN
560      CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,1, &
561        & XCoord, &
562        & Err)
563      CALL CMISSField_ParameterSetGetNode(GeometricField,CMISS_FIELD_U_VARIABLE_TYPE,CMISS_FIELD_VALUES_SET_TYPE,1,1,NODE,2, &
564        & YCoord, &
565        & Err)
566      IF(ABS(XCoord)<1.0E-6_CMISSDP) THEN
567        CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,1, &
568          & CMISS_BOUNDARY_CONDITION_FIXED,XCoord,Err)
569        WRITE(*,*) "FIXING NODE",NODE,"IN X DIRECTION"
570        X_FIXED=.TRUE.
571      ENDIF
572      IF(ABS(YCoord)<1.0E-6_CMISSDP) THEN
573        CALL CMISSBoundaryConditions_SetNode(BoundaryConditions,DependentField,CMISS_FIELD_U_VARIABLE_TYPE,1,1,NODE,2, &
574          & CMISS_BOUNDARY_CONDITION_FIXED,YCoord,Err)
575        WRITE(*,*) "FIXING NODE",NODE,"IN Y DIRECTION"
576        Y_FIXED=.TRUE.
577      ENDIF
578    ENDIF
579  ENDDO
580  CALL MPI_REDUCE(X_FIXED,X_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
581  CALL MPI_REDUCE(Y_FIXED,Y_OKAY,1,MPI_LOGICAL,MPI_LOR,0,MPI_COMM_WORLD,MPI_IERROR)
582  IF(ComputationalNodeNumber==0) THEN
583    IF(.NOT.(X_OKAY.AND.Y_OKAY)) THEN
584      WRITE(*,*) "Free body motion could not be prevented!"
585      CALL CMISSFinalise(Err)
586      STOP
587    ENDIF
588  ENDIF
589
590  ! END ASSIGN BOUNDARY CONDITIONS ===============
591
592  CALL CMISSSolverEquations_BoundaryConditionsCreateFinish(SolverEquations,Err)
593
594  !Solve for many timesteps while incrementing inner pressure
595  DO time_step=1,3
596    !Solve problem
597    CALL CMISSProblem_Solve(Problem,Err)
598
599    !Output solution
600    CALL CMISSFields_Initialise(Fields,Err)
601    CALL CMISSFields_Create(Region,Fields,Err)
602    WRITE(fileNumber,'(I3.3)') time_step
603    fileName='./outputs/QuadraticEllipsoid'//trim(fileNumber)
604    CALL CMISSFields_NodesExport(Fields,trim(fileName),"FORTRAN",Err)
605    CALL CMISSFields_ElementsExport(Fields,trim(fileName),"FORTRAN",Err)
606    CALL CMISSFields_Finalise(Fields,Err)
607
608    !Apply inner surface pressure
609    !NOTE: Surface pressure goes into pressure_values_set_type of the DELUDELN type
610    DO NN=1,SIZE(InnerSurfaceNodes,1)
611      NODE=InnerSurfaceNodes(NN)
612      CALL CMISSDecomposition_NodeDomainGet(Decomposition,NODE,1,NodeDomain,Err)
613      IF(NodeDomain==ComputationalNodeNumber) THEN
614        !Set the field directly
615        CALL CMISSField_ParameterSetUpdateNode(DependentField,CMISS_FIELD_DELUDELN_VARIABLE_TYPE, &
616          & CMISS_FIELD_PRESSURE_VALUES_SET_TYPE,1,1,&
617          & NODE, ABS(InnerNormalXi),INNER_PRESSURE*REAL(time_step+1,CMISSDP),Err)
618      ENDIF
619    ENDDO
620  ENDDO
621
622  CALL CMISSFinalise(Err)
623
624  WRITE(*,'(A)') "Program successfully completed."
625
626  STOP
627
628END PROGRAM QUADRATICELLIPSOIDCOSTAEXAMPLE
629