diff --git a/files/MonteCarloTest/infile_sphere.txt b/files/MonteCarloTest/infile_sphere.txt deleted file mode 100644 index fdbfb3019..000000000 --- a/files/MonteCarloTest/infile_sphere.txt +++ /dev/null @@ -1,42 +0,0 @@ -sphere /* data file */ -2 /* tissptr->num_layers */ -1.000 /* photptr->layerprops[0].n idx of outside medium */ -1.4 /* photptr->layerprops[1].n idx of layer 1 */ -26.75 /* photptr->layerprops[1].mus */ -0.088 /* photptr->layerprops[1].mua */ -0.80 /* photptr->layerprops[1].g */ -10.0 /* photptr->layerprops[1].d */ -1.4 /* photptr->layerprops[1].n idx of layer 1 */ -214 /* photptr->layerprops[1].mus */ -0.088 /* photptr->layerprops[1].mua */ -0.80 /* photptr->layerprops[2].g */ -0.0 /* photptr->layerprops[2].d */ -1.000 /* photptr->layerprops[3].n idx of out bot med */ -r /* beam type f=uniform circular, g=gaussian circular, r=rectangle -0.0000000 /* photptr->offset (>0 or<0) or the beam along x axis*/ -0.05 /* photptr->beam_radius or max x of the rectangle*/ -1.38 /* photptr->src_NA */ -10 /* tissptr->nr */ -0.100 /* tissptr->dr */ */ -10 /* tissptr->nz */ -0.100 /* tissptr->dz */ -80 /* tissptr->nx */ //quadrato 8cm*8cm con passo 1/2mm -0.05 /* tissptr->dx */ */ -80 /* tissptr->ny */ -0.05 /* tissptr->dy */ -100000 /* tissptr->num_photons */ -10 /* tissptr->nt */ -10 /* tissptr->dt */ -3 /* pertptr->do_ellip_layer 0=no pert, 1=ellip, 2=layer, 3=ellip no pMC */ -0.0 /* pertptr->ellip_x */ -0.0 /* pertptr->ellip_y */ -0.48 /* pertptr->ellip_z */ -0.23 /* pertptr->ellip_rad_x */ -0.23 /* pertptr->ellip_rad_y */ -0.23 /* pertptr->ellip_rad_z */ -0.0 /* pertptr->layer_z_min THESE NEED TO MATCH PERT REGION */ -0.02 /* pertptr->layer_z_max */ -1 /* pertptr->nr */ -1 /* pertptr->reflect_flag */ -0.0 /* pertptr->det_ctr[0] */ -3.0 /* pertptr->det_rad */ diff --git a/files/MonteCarloTest/infile_sphere.xml b/files/MonteCarloTest/infile_sphere.xml deleted file mode 100644 index 5e38144de..000000000 --- a/files/MonteCarloTest/infile_sphere.xml +++ /dev/null @@ -1,75 +0,0 @@ - - - sphere - 0 - 3 - 3 - 0.23 - 0.23 - 0.23 - 0 - 0 - 0.48 - 0.02 - 0 - 1 - true - 0.0 - 0.05 - r - - - 0 - 0 - 0 - 0 - 0 - 1 - 0 - 0 - - - 0 - 10 - 0.8 - 0.088 - 26.75 - 1.4 - 0 - 0 - - - 0 - 0 - 0.8 - 0.088 - 214 - 1.4 - 0 - 0 - - - 0 - 0 - 0 - 0 - 0 - 1.4 - 0 - 0 - - - 1.38 - 0.1 - 10 - 0.05 - 0.05 - 0.1 - 10 - 10 - 2 - 10000 - 80 - 80 - 10 - \ No newline at end of file diff --git a/src/.nuget/NuGet.Config b/src/.nuget/NuGet.Config deleted file mode 100644 index 79e7784da..000000000 --- a/src/.nuget/NuGet.Config +++ /dev/null @@ -1,6 +0,0 @@ - - - - - - \ No newline at end of file diff --git a/src/Vts.FemModeling/MGRTE/2D/BoundaryCoupling.cs b/src/Vts.FemModeling/MGRTE/2D/BoundaryCoupling.cs deleted file mode 100644 index 35f873db4..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/BoundaryCoupling.cs +++ /dev/null @@ -1,48 +0,0 @@ -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Boundary reflection or refraction coupling structure - /// - public struct BoundaryCoupling - { - /// - /// reflection contribution to incoming flux from outgoing flux: Ri[ne][ns][2] --- interpolated direction - /// - public int[][][] Ri; - - /// - /// reflection contribution to outgoing flux from boundary source: Ro[ne][ns][2] --- interpolated weight - /// - public int[][][] Ro; - - /// - /// refraction contribution to incoming flux from boundary source: Si[ne][ns][2] --- interpolated direction - /// - public int[][][] Si; - - /// - /// refraction contribution to outgoing flux from outgoing flux: So[ne][ns][2] --- interpolated weight - /// - public int[][][] So; - - /// - /// reflection contribution to incoming flux from outgoing flux: Ri2[ne][ns][2] --- interpolated weight - /// - public double[][][] Ri2; - - /// - /// reflection contribution to outgoing flux from boundary source: Ro2[ne][ns][2] --- interpolated direction - /// - public double[][][] Ro2; - - /// - /// refraction contribution to incoming flux from boundary source: Si2[ne][ns][2] --- interpolated weight - /// - public double[][][] Si2; - - /// - /// refraction contribution to outgoing flux from outgoing flux: So2[ne][ns][2] --- interpolated direction - /// - public double[][][] So2; - }; -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/IntPointSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/IntPointSourceInput.cs deleted file mode 100644 index 5f282702b..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/IntPointSourceInput.cs +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/Measurement.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/Measurement.cs deleted file mode 100644 index f1ec9f436..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/Measurement.cs +++ /dev/null @@ -1,49 +0,0 @@ -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Measurement structure - /// - public struct Measurement - { - /// - /// fluence at each node - /// - public double[] fluence; - - /// - /// radiance at each node - /// - public double[][] radiance; - - /// - /// rectangular grid array - /// - public double[][] uxy; - - /// - /// x coordinates - /// - public double[] xloc; - - /// - /// z coordinates - /// - public double[] zloc; - - /// - /// Intensity - /// - public double[] inten; - - /// - /// dx spacing between each grid - /// - public double[] dx; - - /// - /// dz spacing between each grid - /// - public double[] dz; - } -} - diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/AngularMesh.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/AngularMesh.cs deleted file mode 100644 index fb5d5d877..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/AngularMesh.cs +++ /dev/null @@ -1,21 +0,0 @@ -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Angular mesh structure - /// - public struct AngularMesh - { - /// - /// number of angular nodes - /// - public int Ns; - /// - /// angular coordinates: a[ns][3] with a[i][3]:=(cos(theta(i)), sin(theta(i)), theta(i)) - /// - public double[][] Ang; - /// - /// angular weights: w[ns][ns] - /// - public double[][][] W; - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/MeshSimulationOptions.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/MeshSimulationOptions.cs deleted file mode 100644 index 19c9edb42..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/MeshSimulationOptions.cs +++ /dev/null @@ -1,118 +0,0 @@ -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// MG-RTE simulation parameters - /// - public class MeshSimulationOptions - { - /// - /// Constructor for MG-RTE simulation parameters - /// - /// Refractive index of the external medium - /// The residual value of the iteration for stopping criterion - /// The choice of multigrid method - /// The maximum number of iteration on the finest level in FMG - /// The number of pre-iterations with the suggested value "3". See paper. - /// The number of post-iterations with the suggested value "3". See paper. - /// The number of multigrid cycles on each level except the finest level in FMG with the - /// suggested value "1".See paper. - /// The indicator of full multigrid method (FMG) with the suggested value "1". See paper - /// Starting SmeshLevel (should be less than smeshLevel) - /// Starting AmeshLevel (should be less than ameshLevel) - public MeshSimulationOptions( - double nExternal, - double convTolerance, - int methodMg, - int nIterations, - int nPreIterations, - int nPostIterations, - int nCycle, - int fullMg, - int startingSmeshLevel, - int startingAmeshLevel - ) - { - NExternal = nExternal; - ConvTolerance = convTolerance; - MethodMg = methodMg; - NIterations = nIterations; - NPreIterations = nPreIterations; - NPostIterations = nPostIterations; - NCycle = nCycle; - FullMg = fullMg; - StartingSmeshLevel = startingSmeshLevel; - StartingAmeshLevel = startingAmeshLevel; - } - - /// - /// Initializes a new instance of MG-RTE simulation parameters - /// - /// Refractive index of the external medium - /// The residual value of the iteration for stopping criterion - /// The choice of multigrid method - /// The maximum number of iteration on the finest level in FMG - public MeshSimulationOptions( - double nExternal, - double convTolerance, - int methodMg, - int nIterations) - : this(nExternal, convTolerance, methodMg, nIterations, 3, 3, 1, 1, 1, 1 ) { } - - /// - /// Default constructor for MG-RTE simulation parameters - /// - public MeshSimulationOptions() - : this(1.0, 1.0e-4, 6, 100, 3, 3, 1, 1, 1, 1) { } - - /// - /// Refractive index of the external medium - /// - public double NExternal { get; set; } - - /// - /// The residual value of the iteration for stopping criterion - /// - public double ConvTolerance { get; set; } - - /// - /// The choice of multigrid method - /// - public int MethodMg { get; set; } - - /// - /// The maximum number of iteration on the finest level in FMG - /// - public int NIterations { get; set; } - - /// - /// The number of pre-iterations with the suggested value "3", see the paper. - /// - public int NPreIterations { get; set; } - - /// - /// The number of post-iterations with the suggested value "3", see the paper. - /// - public int NPostIterations { get; set; } - - /// - /// The number of multigrid cycles on each level except the finest level in - /// FMG with the suggested value "1", see the paper. - /// - public int NCycle { get; set; } - - /// - /// The indicator of full multigrid method (FMG) with the suggested value "1" - /// - public int FullMg { get; set; } - - /// - /// Starting SmeshLevel (should be less than smeshLevel) - /// - public int StartingSmeshLevel { get; set; } - - /// - /// Starting AmeshLevel (should be less than ameshLevel) - /// - public int StartingAmeshLevel { get; set; } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/SpatialMesh.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/SpatialMesh.cs deleted file mode 100644 index 70031f4d2..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/SpatialMesh.cs +++ /dev/null @@ -1,113 +0,0 @@ -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Spatial mesh structure - /// - public struct SpatialMesh - { - /// - /// number of node points - /// - public int Np; - - /// - /// number of boundary edges - /// - public int Ne; - - /// - /// number of triangles - /// - public int Nt; - - /// - /// tissue region Region[np] - /// - public int []Region; - - /// - /// nodal coordinates: p[np][2] - /// - public double[][] P; - - /// - /// nodes contained in one triangle: t[nt][3] - /// - public int[][] T; - - /// - /// nodes contained in one boundary edge: e[ne][4] - /// - public int[][] E; - - /// - /// sweep ordering: so[ns][nt] - /// - public int[][] So; - - /// - /// center of triangles: c[nt][2] (see "initialization") - /// - public double[][] C; - - /// - /// center of edges: ec[ne][2] (see "initialization") - /// - public double[][] Ec; - - /// - /// area of triangles: a[nt] (see "initialization") - /// - public double[] A; - - /// - /// triangles adjacent to one node: p2[np][p2[np][0]+1] (see "initialization") - /// - public int[][] P2; - - /// - /// spatial position mapping between coarse and fine mesh: smap[nt_c][smap[nt_c][0]+1] (see "spatialmapping") - /// - public int[][] Smap; - - /// - /// spatial coarse-to-fine nodal-value mapping: cf[nt_c][3][smap[nt_c][0]][3] (see "spatialmapping2") - /// - public double[][][][] Cf; - - /// - /// spatial fine-to-coarse nodal-value mapping: fc[nt_c][3][smap[nt_c][0]][3] (see "spatialmapping2") - /// - public double[][][][] Fc; - - /// - /// boundary information: ori[ne] (see "boundary") - /// - public int[] Ori; - - /// - /// boundary information: e2[ne][2] (see "boundary") - /// - public int[][] E2; - - /// - /// boundary information: so2[nt][3] (see "boundary") - /// - public int[][] So2; - - /// - /// boundary information: n[ne][2] (see "boundary") - /// - public double[][] N; - - /// - /// edge upwind or outgoing flux: bd[ns][nt][3*3] (see "edgeterm") - /// - public int[][][] Bd; - - /// - /// edge upwind or outgoing flux: bd[ns][nt][3] (see "edgeterm") - /// - public double[][][] Bd2; - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/SquareMeshDataInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/SquareMeshDataInput.cs deleted file mode 100644 index 57567a4d9..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/MeshInputs/SquareMeshDataInput.cs +++ /dev/null @@ -1,42 +0,0 @@ -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Angular and spatial mesh data for a square mesh - /// - public class SquareMeshDataInput - { - /// - /// Constructor for mesh input data - /// - /// The finest layer of angular mesh generation - /// The finest layer of spatial mesh generation - /// Length of the square mesh - public SquareMeshDataInput(int aMeshLevel, int sMeshLevel, double sideLength) - { - AMeshLevel = aMeshLevel; - SMeshLevel = sMeshLevel; - SideLength = sideLength; - } - - /// - /// Default constructor for mesh input data - /// - public SquareMeshDataInput() - : this(5, 3, 10.0) { } - - /// - /// The finest layer of angular mesh generation - /// - public int AMeshLevel { get; set; } - - /// - /// The finest layer of spatial mesh generation - /// - public int SMeshLevel { get; set; } - - /// - /// Length of the square mesh - /// - public double SideLength { get; set; } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SimulationInputs.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SimulationInputs.cs deleted file mode 100644 index 3b728beca..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SimulationInputs.cs +++ /dev/null @@ -1,66 +0,0 @@ -using Vts.FemModeling.MGRTE._2D.SourceInputs; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Tissues; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// FEM 2D MG_RTE Simulation inputs - /// - public class SimulationInput - { - /// - /// Input data for spatial and angular mesh - /// - public SquareMeshDataInput MeshDataInput; - /// - /// Mesh simulation parameters - /// - public MeshSimulationOptions SimulationOptionsInput; - /// - /// Specifying external source - /// - public IExtFemSourceInput ExtSourceInput; - /// - /// Specifying internal source - /// - public IIntFemSourceInput IntSourceInput; - /// - /// Specifying tissue definition - /// - public ITissueInput TissueInput; - - /// - /// General constructor for simulation inputs - /// - /// Input data for spatial and angular mesh - /// Mesh simulation options input - /// Specifying external source - /// Specifying internal source - /// Specifying tissue definition - public SimulationInput( - SquareMeshDataInput meshDataInput, - MeshSimulationOptions simulationOptionsInput, - IExtFemSourceInput extSourceInput, - IIntFemSourceInput intSourceInput, - ITissueInput tissueInput - ) - { - MeshDataInput = meshDataInput; - SimulationOptionsInput = simulationOptionsInput; - ExtSourceInput = extSourceInput; - IntSourceInput = intSourceInput; - TissueInput = tissueInput; - } - - /// - /// Default constructor for simulation inputs - /// - public SimulationInput() - : this(new SquareMeshDataInput(), - new MeshSimulationOptions(), - new ExtPointSourceInput(), - null, - new MultiEllipsoidTissueInput()) { } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/ExternalSourceInputs/ExtLineSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/ExternalSourceInputs/ExtLineSourceInput.cs deleted file mode 100644 index b99b818e1..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/ExternalSourceInputs/ExtLineSourceInput.cs +++ /dev/null @@ -1,55 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.SourceInputs -{ - /// - /// External line source - /// - public class ExtLineSourceInput : IExtFemSourceInput - { - /// - /// General constructor for external line source - /// - /// starting point (x,z) - /// end point (x,z) - /// Theta Range - public ExtLineSourceInput( - DoubleRange start, - DoubleRange end, - DoubleRange thetaRange) - { - SourceType = FemSourceType.ExtPointSource; - Start = start; - End = end; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for external line source - /// - public ExtLineSourceInput() - : this( - new DoubleRange(0.0, 0.0), - new DoubleRange(0.0, 0.0), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Ext source type - /// - public FemSourceType SourceType { get; set; } - /// - /// Starting cooordinates (x,z) - /// - public DoubleRange Start { get; set; } - /// - /// Ending cooordinates (x,z) - /// - public DoubleRange End { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/ExternalSourceInputs/ExtPointSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/ExternalSourceInputs/ExtPointSourceInput.cs deleted file mode 100644 index 07a930961..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/ExternalSourceInputs/ExtPointSourceInput.cs +++ /dev/null @@ -1,47 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.SourceInputs -{ - /// - /// External point source - /// - public class ExtPointSourceInput : IExtFemSourceInput - { - /// - /// General constructor for external point source - /// - /// Launch point (x,z) - /// Theta Range - public ExtPointSourceInput( - DoubleRange launchPoint, - DoubleRange thetaRange) - { - SourceType = FemSourceType.ExtPointSource; - LaunchPoint = launchPoint; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for external point source - /// - public ExtPointSourceInput() - : this( - new DoubleRange(0.0, 0.0), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Ext source type - /// - public FemSourceType SourceType { get; set; } - /// - /// Launching cooordinates (x,z) - /// - public DoubleRange LaunchPoint { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DCircularSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DCircularSourceInput.cs deleted file mode 100644 index beb685f18..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DCircularSourceInput.cs +++ /dev/null @@ -1,55 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.SourceInputs -{ - /// - /// Internal 2D Circular source - /// - public class Int2DCircularSourceInput : IIntFemSourceInput - { - /// - /// General constructor for 2D circular source - /// - /// Radius - /// Center (x,z) - /// Theta Range - public Int2DCircularSourceInput( - double radius, - DoubleRange center, - DoubleRange thetaRange) - { - SourceType = FemSourceType.Int2DCircularSource; - Radius = radius; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D circular source - /// - public Int2DCircularSourceInput() - : this( - 0.25, - new DoubleRange(0.0, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Internal source type - /// - public FemSourceType SourceType { get; set; } - /// - /// Radius of the circular source - /// - public double Radius { get; set; } - /// - /// Center cooordinates (x,z) - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DEllipticalSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DEllipticalSourceInput.cs deleted file mode 100644 index f4b4d77d3..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DEllipticalSourceInput.cs +++ /dev/null @@ -1,63 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.SourceInputs -{ - /// - /// Internal 2D Elliptical source - /// - public class Int2DEllipticalSourceInput : IIntFemSourceInput - { - /// - /// General constructor for 2D elliptical source - /// - /// a Parameter - /// b parameter - /// Center (x,z) - /// Theta Range - public Int2DEllipticalSourceInput( - double aParameter, - double bParameter, - DoubleRange center, - DoubleRange thetaRange) - { - SourceType = FemSourceType.Int2DEllipticalSource; - AParameter = aParameter; - BParameter = bParameter; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D elliptical source - /// - public Int2DEllipticalSourceInput() - : this( - 0.25, - 0.125, - new DoubleRange(0.0, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Internal source type - /// - public FemSourceType SourceType { get; set; } - /// - /// a Parameter of the Elliptical source - /// - public double AParameter { get; set; } - /// - /// b Parameter of the Elliptical source - /// - public double BParameter { get; set; } - /// - /// Center cooordinates (x,z) - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DPointSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DPointSourceInput.cs deleted file mode 100644 index 7e629bce8..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DPointSourceInput.cs +++ /dev/null @@ -1,47 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.SourceInputs -{ - /// - /// Internal 2D Point source - /// - public class Int2DPointSourceInput : IIntFemSourceInput - { - /// - /// General constructor for 2D point source - /// - /// Center (x,z) - /// Theta Range - public Int2DPointSourceInput( - DoubleRange center, - DoubleRange thetaRange) - { - SourceType = FemSourceType.Int2DPointSource; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D point source - /// - public Int2DPointSourceInput() - : this( - new DoubleRange(0, 5.0), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Internal source type - /// - public FemSourceType SourceType { get; set; } - /// - /// Center cooordinates (x,z) - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DRectangularSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DRectangularSourceInput.cs deleted file mode 100644 index 73dafb641..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/DataStructures/SourceInputs/InternalSourceInputs/Int2DRectangularSourceInput.cs +++ /dev/null @@ -1,63 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.SourceInputs -{ - /// - /// Internal 2D Rectangular source - /// - public class Int2DRectangularSourceInput : IIntFemSourceInput - { - /// - /// General constructor for 2D rectangular source - /// - /// X Length - /// Z Length - /// Center (x,z) - /// Theta Range - public Int2DRectangularSourceInput( - double xLength, - double zHeight, - DoubleRange center, - DoubleRange thetaRange) - { - SourceType = FemSourceType.Int2DRectangularSource; - XLength = xLength; - ZHeight = zHeight; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D rectangular source - /// - public Int2DRectangularSourceInput() - : this( - 0.25, - 0.125, - new DoubleRange(0.0, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// External source type - /// - public FemSourceType SourceType { get; set; } - /// - /// Length of the rectangular source - /// - public double XLength { get; set; } - /// - /// height of the rectangular source - /// - public double ZHeight { get; set; } - /// - /// Center cooordinates (x,z) - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Enums.cs b/src/Vts.FemModeling/MGRTE/2D/Enums.cs deleted file mode 100644 index 1a8cad151..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Enums.cs +++ /dev/null @@ -1,102 +0,0 @@ - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Multigrid Type - /// - public enum WhichMultiGridType - { - /// - /// Angular MG - /// - MgA = 1, - /// - /// Spatial MG - /// - MgS, - /// - /// MG Method 1 - /// - Mg1, - /// - /// MG Method 2 - /// - Mg2, - /// - /// MG Method 3 - /// - Mg3, - /// - /// Angular -> Spatial MG - /// - Mg4A, - /// - /// Spatial -> Angular MG - /// - Mg4S - } - - /// - /// Source type - /// - public enum FemSourceType - { - /// - /// External line Source - /// - ExtLineSource, - /// - /// External Point source - /// - ExtPointSource, - /// - /// Internal 2D Circular Source - /// - Int2DCircularSource, - /// - /// Internal 2D Elliptical Source - /// - Int2DEllipticalSource, - /// - /// Internal 2D Rectangular Source - /// - Int2DRectangularSource, - /// - /// Internal 2D Point Source - /// - Int2DPointSource - } - - /// - /// Tally type - /// - public enum TallyType - { - /// - /// Internal fluence distribution - /// - FluenceOfXAndZ, - } - - ///// - ///// Tissue type - ///// - //public enum TissueType - //{ - // /// - // /// Multilayer tissue type. Includes homogeneous tissues. - // /// - // MultiLayer, - //} - - /// - /// Inclusion type - /// - public enum InclusionType - { - /// - /// Circular inclusion in multi layer - /// - MultiLayerCircular, - } -} \ No newline at end of file diff --git a/src/Vts.FemModeling/MGRTE/2D/Factory/FemSourceFactory.cs b/src/Vts.FemModeling/MGRTE/2D/Factory/FemSourceFactory.cs deleted file mode 100644 index bb50cbfbe..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Factory/FemSourceFactory.cs +++ /dev/null @@ -1,84 +0,0 @@ -using System; -using Vts.FemModeling.MGRTE._2D.DataStructures; -using Vts.FemModeling.MGRTE._2D.SourceInputs; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Instantiates appropriate source - /// - public static class FemSourceFactory - { - /// - /// Get External source - /// - /// External FEM source - /// selected source - public static IExtSource GetExtSource(IExtFemSourceInput input) - { - switch (input.SourceType) - { - case FemSourceType.ExtPointSource: - var epsInput = (ExtPointSourceInput)input; - return new ExtPointSource( - epsInput.LaunchPoint, - epsInput.ThetaRange); - - case FemSourceType.ExtLineSource: - var elsInput = (ExtLineSourceInput)input; - return new ExtLineSource( - elsInput.Start, - elsInput.End, - elsInput.ThetaRange); - - default: - throw new NotImplementedException( - "Problem generating IExtSource instance. Check that IExtFemSourceInput has a matching IExtSource definition."); - } - } - - /// - /// Get Internal source - /// - /// Internal FEM source - /// selected source - public static IIntSource GetIntSource(IIntFemSourceInput input) - { - switch (input.SourceType) - { - case FemSourceType.Int2DPointSource: - var ipsInput = (Int2DPointSourceInput)input; - return new Int2DPointSource( - ipsInput.Center, - ipsInput.ThetaRange); - - case FemSourceType.Int2DCircularSource: - var icsInput = (Int2DCircularSourceInput)input; - return new Int2DCircularSource( - icsInput.Radius, - icsInput.Center, - icsInput.ThetaRange); - - case FemSourceType.Int2DEllipticalSource: - var iesInput = (Int2DEllipticalSourceInput)input; - return new Int2DEllipticalSource( - iesInput.AParameter, - iesInput.BParameter, - iesInput.Center, - iesInput.ThetaRange); - - case FemSourceType.Int2DRectangularSource: - var irsInput = (Int2DRectangularSourceInput)input; - return new Int2DRectangularSource( - irsInput.XLength, - irsInput.ZHeight, - irsInput.Center, - irsInput.ThetaRange); - - default: - throw new NotImplementedException( - "Problem generating IIntSource instance. Check that IIntFemSourceInput has a matching IIntSource definition."); - } - } - } -} \ No newline at end of file diff --git a/src/Vts.FemModeling/MGRTE/2D/Initialization.cs b/src/Vts.FemModeling/MGRTE/2D/Initialization.cs deleted file mode 100644 index c4902ff3f..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Initialization.cs +++ /dev/null @@ -1,582 +0,0 @@ -using System; -using Vts.FemModeling.MGRTE._2D.DataStructures; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Initialization class - /// - public class Initialization - { - public static void Initial( - ref AngularMesh[] amesh, - ref SpatialMesh[] smesh, - ref double[][][][] flux, - ref double[][][][] d, - ref double[][][][] RHS, - ref double[][][][] q, - ref int[][] noflevel, - ref BoundaryCoupling[] b, - int level, - int mgMethod, - int vacuum, - double nTissue, - double nExt, - int aMeshLevel, - int aMeshLevel0, - int sMeshLevel, - int sMeshLevel0, - double[][][] ua, - double[][][] us, - MultiGridCycle Mgrid) - { - - - //Purpose: in this function, we load the following four ".txt" files: - // 1) "amesh.txt": "ns", "a" and "w" for angular mesh - // 3) "ua.txt": absorption coefficient - // 4) "us.txt": scattering coefficient - - // we first load the mesh files: - // 1.1. "amesh.txt": "ns", "a" and "w" for angular mesh - - // and then we compute the following for spatial mesh: - // 2.1. "smap", "cf" and "fc" - - // and finally malloc the following for multigrid algorithm: - // 3.1. "noflevel" - // 3.2. "ua", "us", "flux", "RHS", "d" and "q" - // 3.3. "b" - - int tempsize = 20, tempsize2 = 20; - int i, j, k, m, n, nt, np, ne, ns, da = aMeshLevel - aMeshLevel0, ds = sMeshLevel - sMeshLevel0; - double x1, x2, x3, y1, y2, y3; - int[][] t; - int[][] e; - int[][] p2; - int[][] smap; - double[][] p; - - - // 2.1. compute "c", "ec", "a" and "p2" - // p2[np][p2[np][0]+1]: triangles adjacent to one node - // For the 2nd index, the 1st element is the total number of triangles adjacent to this node, - // the corresponding triangles are saved from the 2nd. - // Example: assume triangles [5 16 28 67] are adjacent to the node 10, then - // p2[10][0]=4, p2[10][1]=5, p2[10][2]=16, p2[10][3]=28 and p2[10][4]=67. - for (i = 0; i <= sMeshLevel; i++) - { - p = smesh[i].P; t = smesh[i].T; nt = smesh[i].Nt; np = smesh[i].Np; e = smesh[i].E; ne = smesh[i].Ne; - - smesh[i].C = new double[nt][]; - for (j = 0; j < nt; j++) - { - smesh[i].C[j] = new double[2]; - for (k = 0; k < 2; k++) - { - smesh[i].C[j][k] = (p[t[j][0]][k] + p[t[j][1]][k] + p[t[j][2]][k]) / 3; - }// center of triangle - } - - smesh[i].Ec = new double[ne][]; - for (j = 0; j < ne; j++) - { - smesh[i].Ec[j] = new double[2]; ; - for (k = 0; k < 2; k++) - { smesh[i].Ec[j][k] = (p[e[j][1]][k] + p[e[j][2]][k]) / 2; }// center of edge - } - - smesh[i].A = new double[nt]; ; - for (j = 0; j < nt; j++) - { - x1 = p[t[j][0]][0]; y1 = p[t[j][0]][1]; - x2 = p[t[j][1]][0]; y2 = p[t[j][1]][1]; - x3 = p[t[j][2]][0]; y3 = p[t[j][2]][1]; - smesh[i].A[j] = MathFunctions.Area(x1, y1, x2, y2, x3, y3);//area of triangle - } - - p2 = new int[np][]; - for (j = 0; j < np; j++) - { - p2[j] = new int[tempsize]; - // tempsize is the initial length of the second index of "p2", and it may cause problem if it is too small. - p2[j][0] = 0; - for (k = 1; k < tempsize; k++) - { - p2[j][k] = -1; - } - } - - for (j = 0; j < nt; j++) - { - for (k = 0; k < 3; k++) - { - p2[t[j][k]][0] += 1; - p2[t[j][k]][p2[t[j][k]][0]] = j; - } - } - - for (j = 0; j < np; j++) - { - if (p2[j][0] + 1 > tempsize) - { - Console.WriteLine("WARNING: tempsize for p2 is too small!!\n"); - goto stop; - } - } - - smesh[i].P2 = new int[np][]; - for (j = 0; j < np; j++) - { - smesh[i].P2[j] = new int[(p2[j][0] + 1)]; - for (k = 0; k <= p2[j][0]; k++) - { smesh[i].P2[j][k] = p2[j][k]; } - } - } - - // 2.2. compute "smap", "cf" and "fc" - // For the data structure of "smap", see "spatialmapping"; - // For the data structure of "cf" and "fc", see "spatialmapping2". - // Note: those arrays are not defined on the coarsest spatial mesh (i=0), - // since they are always saved on the fine level instead of the coarse level. - for (i = 1; i <= sMeshLevel; i++) - { - smap = new int[smesh[i - 1].Nt][]; - for (j = 0; j < smesh[i - 1].Nt; j++) - { - smap[j] = new int[tempsize2]; - // tempsize2 is the initial length of the second index of "smap", and it may cause problem if it is too small. - smap[j][0] = 0; - for (k = 1; k < tempsize2; k++) - { smap[j][k] = -1; } - } - Mgrid.SpatialMapping(smesh[i - 1], smesh[i], smap); - - for (j = 0; j < smesh[i - 1].Nt; j++) - { - if (smap[j][0] > tempsize2 - 1) - { - Console.WriteLine("WARNING: tempsize2 for smap is too small!!\n"); - goto stop; - } - } - - smesh[i].Smap = new int[smesh[i - 1].Nt][]; - for (j = 0; j < smesh[i - 1].Nt; j++) - { - smesh[i].Smap[j] = new int[smap[j][0] + 1]; - for (k = 0; k <= smap[j][0]; k++) - { smesh[i].Smap[j][k] = smap[j][k]; } - } - - - smesh[i].Cf = new double[smesh[i - 1].Nt][][][]; - for (j = 0; j < smesh[i - 1].Nt; j++) - { - smesh[i].Cf[j] = new double[3][][]; - for (k = 0; k < 3; k++) - { - smesh[i].Cf[j][k] = new double[smesh[i].Smap[j][0]][]; - for (m = 0; m < smesh[i].Smap[j][0]; m++) - { smesh[i].Cf[j][k][m] = new double[3]; ; } - } - } - smesh[i].Fc = new double[smesh[i - 1].Nt][][][]; - for (j = 0; j < smesh[i - 1].Nt; j++) - { - smesh[i].Fc[j] = new double[3][][]; - for (k = 0; k < 3; k++) - { - smesh[i].Fc[j][k] = new double[smesh[i].Smap[j][0]][]; - for (m = 0; m < smesh[i].Smap[j][0]; m++) - { smesh[i].Fc[j][k][m] = new double[3]; ; } - } - } - Mgrid.SpatialMapping2(smesh[i - 1], smesh[i], smesh[i].Smap, smesh[i].Cf, smesh[i].Fc, tempsize2); - } - - // 2.3. compute "e", "e2", "so2", "n" and "ori" - // For the data structure of "eo", "e2","so2", "n" and "ori", see "boundary". - for (i = 0; i <= sMeshLevel; i++) - { - smesh[i].E2 = new int[smesh[i].Ne][]; - for (j = 0; j < smesh[i].Ne; j++) - { smesh[i].E2[j] = new int[2]; } - smesh[i].So2 = new int[smesh[i].Nt][]; - for (j = 0; j < smesh[i].Nt; j++) - { smesh[i].So2[j] = new int[3]; } - smesh[i].N = new double[smesh[i].Ne][]; - for (j = 0; j < smesh[i].Ne; j++) - { smesh[i].N[j] = new double[2]; } - smesh[i].Ori = new int[smesh[i].Ne]; - - Mgrid.Boundary(smesh[i].Ne, smesh[i].Nt, smesh[i].T, smesh[i].P2, smesh[i].P, smesh[i].E, smesh[i].E2, smesh[i].So2, smesh[i].N, smesh[i].Ori); - } - - // 2.4. compute "bd" and "bd2" - // For the data structure of "bd" and "bd2", see "edgeterm". - for (i = 0; i <= sMeshLevel; i++) - { - smesh[i].Bd = new int[amesh[aMeshLevel].Ns][][]; - for (j = 0; j < amesh[aMeshLevel].Ns; j++) - { - smesh[i].Bd[j] = new int[smesh[i].Nt][]; ; - for (k = 0; k < smesh[i].Nt; k++) - { - smesh[i].Bd[j][k] = new int[9]; - for (m = 0; m < 9; m++) - { smesh[i].Bd[j][k][m] = -1; } - } - } - } - for (i = 0; i <= sMeshLevel; i++) - { - smesh[i].Bd2 = new double[amesh[aMeshLevel].Ns][][]; - for (j = 0; j < amesh[aMeshLevel].Ns; j++) - { - smesh[i].Bd2[j] = new double[smesh[i].Nt][]; - for (k = 0; k < smesh[i].Nt; k++) - { smesh[i].Bd2[j][k] = new double[3]; } - } - } - for (i = 0; i <= sMeshLevel; i++) - { - for (j = 0; j < amesh[aMeshLevel].Ns; j++) - { - Mgrid.EdgeTri(smesh[i].Nt, amesh[aMeshLevel].Ang[j], smesh[i].P, smesh[i].P2, smesh[i].T, smesh[i].Bd[j], smesh[i].Bd2[j], smesh[i].So2); - } - } - - // 3.1. compute "noflevel" - // level: the mesh layers for multgrid; the bigger value represents finer mesh. - // noflevel[level][2]: the corresponding angular mesh level and spatial mesh level for each multigrid mesh level - // Example: assume noflevel[i][0]=3 and noflevel[i][1]=2, then - // the spatial mesh level is "3" and the angular mesh level is "2" on the multgrid mesh level "i". - switch (mgMethod) - { - case 1: //AMG - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - for (i = 0; i <= da; i++) - { - noflevel[i][0] = sMeshLevel; - noflevel[i][1] = i + aMeshLevel0; - } - break; - case 2: //SMG - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - for (i = 0; i <= ds; i++) - { - noflevel[i][0] = i + sMeshLevel0; - noflevel[i][1] = aMeshLevel; - } - break; - case 3: //MG1 - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - if (ds > da) - { - for (i = 0; i < ds - da; i++) - { - noflevel[i][0] = i + sMeshLevel0; - noflevel[i][1] = aMeshLevel0; - } - for (i = 0; i <= da; i++) - { - noflevel[i + ds - da][0] = i + ds - da + sMeshLevel0; - noflevel[i + ds - da][1] = i + aMeshLevel0; - } - } - else - { - for (i = 0; i < da - ds; i++) - { - noflevel[i][0] = sMeshLevel0; - noflevel[i][1] = i + aMeshLevel0; - } - for (i = 0; i <= ds; i++) - { - noflevel[i + da - ds][0] = i + sMeshLevel0; - noflevel[i + da - ds][1] = i + da - ds + aMeshLevel0; - } - } - break; - case 4: //MG2 - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - for (i = sMeshLevel0; i <= sMeshLevel; i++) - { - noflevel[i - sMeshLevel0][0] = i; - noflevel[i - sMeshLevel0][1] = aMeshLevel0; - } - for (i = aMeshLevel0 + 1; i <= aMeshLevel; i++) - { - noflevel[i - aMeshLevel0 + sMeshLevel - sMeshLevel0][0] = sMeshLevel; - noflevel[i - aMeshLevel0 + sMeshLevel - sMeshLevel0][1] = i; - } - break; - case 5: //MG3 - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - for (i = aMeshLevel0; i <= aMeshLevel; i++) - { - noflevel[i - aMeshLevel0][0] = sMeshLevel0; - noflevel[i - aMeshLevel0][1] = i; - } - for (i = sMeshLevel0 + 1; i <= sMeshLevel; i++) - { - noflevel[i - sMeshLevel0 + aMeshLevel - aMeshLevel0][0] = i; - noflevel[i - sMeshLevel0 + aMeshLevel - aMeshLevel0][1] = aMeshLevel; - } - break; - case 6: //MG4_a - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - if (ds >= da) - { - for (i = 0; i <= ds - da; i++) - { - noflevel[i][0] = i + sMeshLevel0; - noflevel[i][1] = aMeshLevel0; - } - for (i = 1; i <= da; i++) - { - noflevel[ds - da + 2 * i - 1][0] = sMeshLevel0 + ds - da + i; - noflevel[ds - da + 2 * i - 1][1] = aMeshLevel0 + i - 1; - noflevel[ds - da + 2 * i][0] = sMeshLevel0 + ds - da + i; - noflevel[ds - da + 2 * i][1] = aMeshLevel0 + i; - } - } - else - { - for (i = 0; i <= da - ds; i++) - { - noflevel[i][0] = sMeshLevel0; - noflevel[i][1] = i + aMeshLevel0; - } - for (i = 1; i <= ds; i++) - { - noflevel[da - ds + 2 * i - 1][0] = sMeshLevel0 + i; - noflevel[da - ds + 2 * i - 1][1] = aMeshLevel0 + da - ds + i - 1; - noflevel[da - ds + 2 * i][0] = sMeshLevel0 + i; - noflevel[da - ds + 2 * i][1] = aMeshLevel0 + da - ds + i; - } - } - break; - case 7: //MG4_s - for (i = 0; i <= level; i++) - { noflevel[i] = new int[2]; } - if (ds >= da) - { - for (i = 0; i <= ds - da; i++) - { - noflevel[i][0] = i + sMeshLevel0; - noflevel[i][1] = aMeshLevel0; - } - for (i = 1; i <= da; i++) - { - noflevel[ds - da + 2 * i - 1][0] = sMeshLevel0 + ds - da + i - 1; - noflevel[ds - da + 2 * i - 1][1] = aMeshLevel0 + i; - noflevel[ds - da + 2 * i][0] = sMeshLevel0 + ds - da + i; - noflevel[ds - da + 2 * i][1] = aMeshLevel0 + i; - } - } - else - { - for (i = 0; i <= da - ds; i++) - { - noflevel[i][0] = sMeshLevel0; - noflevel[i][1] = i + aMeshLevel0; - } - for (i = 1; i <= ds; i++) - { - noflevel[da - ds + 2 * i - 1][0] = sMeshLevel0 + i - 1; - noflevel[da - ds + 2 * i - 1][1] = aMeshLevel0 + da - ds + i; - noflevel[da - ds + 2 * i][0] = sMeshLevel0 + i; - noflevel[da - ds + 2 * i][1] = aMeshLevel0 + da - ds + i; - } - } - break; - } - - // 3.2. malloc "ua", "us", "flux", "RHS", "d" and "q" - // ua[level][nt][2]:absorption coefficient - // us[level][nt][2]:scattering coefficient - // flux[level][ns][nt][2]: photon flux - // RHS[level][ns][nt][2]: source term - // d[level][ns][nt][2]: residual - // q[level][2][ns]: boundary source - { - - - for (i = sMeshLevel - 1; i >= 0; i--) - { - ua[i] = new double[smesh[i].Nt][]; - us[i] = new double[smesh[i].Nt][]; - for (j = 0; j < smesh[i].Nt; j++) - { - ua[i][j] = new double[3]; - us[i][j] = new double[3]; - } - Mgrid.FtoC_s2(smesh[i].Nt, ua[i + 1], ua[i], smesh[i + 1].Smap, smesh[i + 1].Fc); - Mgrid.FtoC_s2(smesh[i].Nt, us[i + 1], us[i], smesh[i + 1].Smap, smesh[i + 1].Fc); - } - - - for (n = 0; n <= level; n++) - { - nt = smesh[noflevel[n][0]].Nt; - ns = amesh[noflevel[n][1]].Ns; - RHS[n] = new double[ns][][]; - for (i = 0; i < ns; i++) - { - RHS[n][i] = new double[nt][]; - for (j = 0; j < nt; j++) - { - RHS[n][i][j] = new double[3]; - for (k = 0; k < 3; k++) - { - RHS[n][i][j][k] = 0; - } - } - } - } - - for (n = 0; n <= level; n++) - { - nt = smesh[noflevel[n][0]].Nt; - ns = amesh[noflevel[n][1]].Ns; - d[n] = new double[ns][][]; - for (i = 0; i < ns; i++) - { - d[n][i] = new double[nt][]; - for (j = 0; j < nt; j++) - { - d[n][i][j] = new double[3]; - for (k = 0; k < 3; k++) - { - d[n][i][j][k] = 0; - } - } - } - } - for (n = 0; n <= level; n++) - { - nt = smesh[noflevel[n][0]].Nt; - ns = amesh[noflevel[n][1]].Ns; - flux[n] = new double[ns][][]; - for (i = 0; i < ns; i++) - { - flux[n][i] = new double[nt][]; - for (j = 0; j < nt; j++) - { - flux[n][i][j] = new double[3]; - for (k = 0; k < 3; k++) - { - flux[n][i][j][k] = 0; - } - } - } - } - - for (n = 0; n <= level; n++) - { - ne = smesh[noflevel[n][0]].Ne; - ns = amesh[noflevel[n][1]].Ns; - q[n] = new double[ns][][]; - for (i = 0; i < ns; i++) - { - q[n][i] = new double[ne][]; - for (j = 0; j < ne; j++) - { - q[n][i][j] = new double[2]; - for (k = 0; k < 2; k++) - { q[n][i][j][k] = 0; } - } - } - } - } - - - // 3.3. compute "b" - // For the data structure of "b", see "boundarycoupling". - if (vacuum == 0)// we need "b" only in the presence of refraction index mismatch at the domain boundary. - { - for (i = 0; i <= level; i++) - { - ne = smesh[noflevel[i][0]].Ne; - ns = amesh[noflevel[i][1]].Ns; - - b[i].Ri = new int[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].Ri[j] = new int[ns][]; - for (k = 0; k < ns; k++) - b[i].Ri[j][k] = new int[2]; - } - - b[i].Si = new int[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].Si[j] = new int[ns][]; - for (k = 0; k < ns; k++) - b[i].Si[j][k] = new int[2]; - } - - b[i].Ro = new int[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].Ro[j] = new int[ns][]; - for (k = 0; k < ns; k++) - b[i].Ro[j][k] = new int[2]; - } - - - b[i].So = new int[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].So[j] = new int[ns][]; - for (k = 0; k < ns; k++) - b[i].So[j][k] = new int[2]; - } - - b[i].Ri2 = new double[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].Ri2[j] = new double[ns][]; - for (k = 0; k < ns; k++) - b[i].Ri2[j][k] = new double[2]; - } - - b[i].Ro2 = new double[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].Ro2[j] = new double[ns][]; - for (k = 0; k < ns; k++) - b[i].Ro2[j][k] = new double[2]; - } - - b[i].Si2 = new double[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].Si2[j] = new double[ns][]; - for (k = 0; k < ns; k++) - b[i].Si2[j][k] = new double[2]; - } - - b[i].So2 = new double[ne][][]; - for (j = 0; j < ne; j++) - { - b[i].So2[j] = new double[ns][]; - for (k = 0; k < ns; k++) - b[i].So2[j][k] = new double[2]; - } - - Mgrid.BoundReflection(ns, amesh[noflevel[i][1]].Ang, smesh[noflevel[i][0]], nTissue, nExt, b[i]); - } - } - stop: ; - } - } -} \ No newline at end of file diff --git a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IExtFemSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/Interfaces/IExtFemSourceInput.cs deleted file mode 100644 index aa9b8f13b..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IExtFemSourceInput.cs +++ /dev/null @@ -1,14 +0,0 @@ - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Defines a contract for SourceInput classes. - /// - public interface IExtFemSourceInput - { - /// - /// Type of source - /// - FemSourceType SourceType { get; set; } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IExtSource.cs b/src/Vts.FemModeling/MGRTE/2D/Interfaces/IExtSource.cs deleted file mode 100644 index bc17916d0..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IExtSource.cs +++ /dev/null @@ -1,21 +0,0 @@ -using Vts.FemModeling.MGRTE._2D.DataStructures; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Defines a contract for Ext source classes. - /// - public interface IExtSource - { - /// - /// Assign mesh arrays for an external source - /// - /// - /// - /// - /// - /// - /// - void AssignMeshForExtSource(AngularMesh[] amesh, int ameshLevel, SpatialMesh[] smesh, int smeshLevel, int level, double[][][][] q); - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IInclusionInput.cs b/src/Vts.FemModeling/MGRTE/2D/Interfaces/IInclusionInput.cs deleted file mode 100644 index 99cbe0b10..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IInclusionInput.cs +++ /dev/null @@ -1,20 +0,0 @@ -using Vts.MonteCarlo; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Defines a contract for TissueInput. - /// - public interface IInclusionInput - { - /// - /// Type of tissue - /// - InclusionType InclusionType { get; } - - /// - /// List of tissue regions comprising tissue. - /// - ITissueRegion[] Regions { get; } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IIntFemSourceInput.cs b/src/Vts.FemModeling/MGRTE/2D/Interfaces/IIntFemSourceInput.cs deleted file mode 100644 index 98a5d863e..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IIntFemSourceInput.cs +++ /dev/null @@ -1,14 +0,0 @@ - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Defines a contract for SourceInput classes. - /// - public interface IIntFemSourceInput - { - /// - /// Type of source - /// - FemSourceType SourceType { get; set; } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IIntSource.cs b/src/Vts.FemModeling/MGRTE/2D/Interfaces/IIntSource.cs deleted file mode 100644 index cfeb5b894..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Interfaces/IIntSource.cs +++ /dev/null @@ -1,21 +0,0 @@ -using Vts.FemModeling.MGRTE._2D.DataStructures; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Defines a contract for Internal source classes. - /// - public interface IIntSource - { - /// - /// Assign mesh arrays for an internal source - /// - /// - /// - /// - /// - /// - /// - void AssignMeshForIntSource(AngularMesh[] amesh, int ameshLevel, SpatialMesh[] smesh, int smeshLevel, int level, double[][][][] rhs); - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Interfaces/ITissue.cs b/src/Vts.FemModeling/MGRTE/2D/Interfaces/ITissue.cs deleted file mode 100644 index 121076c99..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Interfaces/ITissue.cs +++ /dev/null @@ -1,48 +0,0 @@ -using System.Collections.Generic; -using Vts.Common; -using Vts.MonteCarlo; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Defines a contract for Tissue classes in Monte Carlo simulation. - /// - public interface ITissue - { - /// - /// AbsorptionWeightingType enum specifier indicating Analog, Discrete Absorption weighting, etc. - /// - AbsorptionWeightingType AbsorptionWeightingType { get; } - - /// - /// photon weight threshold, below which turns on Russian Roulette - /// - double RussianRouletteWeightThreshold { get; } - - /// - /// PhaseFunctionType enum specifier indicating Henyey-Greenstein, Bidirectional, etc. - /// - PhaseFunctionType PhaseFunctionType { get; } - - /// - /// A list of ITissueRegions that describes the entire system. - /// - IList Regions { get; } - - /// - /// The scattering lengths for each tissue region. For discrete absorption weighting and - /// analog, this is based on 1/mut, for continuous it is based on 1/mus. - /// - IList RegionScatterLengths { get; } - - /// - /// Method that gives the current region index within Regions list (above) at the current - /// position of the photon. - /// - /// current location of photon - /// - int GetRegionIndex(Position position); - - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/MathFunctions.cs b/src/Vts.FemModeling/MGRTE/2D/MathFunctions.cs deleted file mode 100644 index 0c18c13d7..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/MathFunctions.cs +++ /dev/null @@ -1,1082 +0,0 @@ -using System; -using Vts.Common; -using Vts.FemModeling.MGRTE._2D.DataStructures; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Tissues; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// Common Math functions - /// - public static class MathFunctions - { - /// - /// Calculate area of a traingle - /// - /// x coordinate 1 - /// y coordinate 1 - /// x coordinate 2 - /// y coordinate 2 - /// x coordinate 3 - /// y coordinate 3 - /// area - public static double Area(double x1, double y1, double x2, double y2, double x3, double y3) - // Purpose: Calculate the area of a triangle - { - return 0.5 * Math.Abs((y2 - y1) * (x3 - x1) - (y3 - y1) * (x2 - x1)); - } - - /// - /// find the minimum from the vector d with the size n. - /// - /// size - /// vector - /// - public static int FindMin(int n, double[] d) - { - int i; - int nmin= 0; - double dmin = d[nmin]; - for (i = 1; i < n; i++) - { - if (d[i] < dmin) - { nmin = i; dmin = d[i]; } - } - return nmin; - } - - /// - /// find two linearly intepolated angles "b" and weights "b2" for the angle "theta_m" - /// - /// theta modified - /// - /// - /// - /// - /// - public static void Intepolation_a(double theta_m, double dtheta, int ns, int[] b, double[] b2, double constant) - { - int theta1, theta2; - double w1, w2; - - theta1 = (int)Math.Floor(theta_m / dtheta) + 1; - w2 = (theta_m - (theta1 - 1) * dtheta) / dtheta; - w1 = 1.0 - w2; - if (theta1 == ns) - { theta2 = 1; } - else - { theta2 = theta1 + 1; } - b[0] = theta1 - 1; b[1] = theta2 - 1; - b2[0] = w1 * constant; b2[1] = w2 * constant; - - } - - /// - /// calculate the distance between two points - /// - /// x coordinate 1 - /// y coordinate 1 - /// x coordinate 2 - /// y coordinate 2 - /// - public static double Length(double x1, double y1, double x2, double y2) - { - return Math.Sqrt((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1)); - } - - /// - /// This function provides angular weights and coordinates on 2D unit circular phase space. - /// - /// Weight - /// Angular array - /// number of angles - /// anisotropy factor - /// - public static void Weight_2D(double[][][] w, double[][] theta, int nAngle, double g, int region) - { - int i, j, n, nj; - double dx, norm; - double[] f; - double[] x; - - n = (nAngle + 2) / 2; - dx = 1 / (double)(n - 1) * 3.14159265; - - f = new double[n]; - x = new double[n]; - - for (i = 0; i < n; i++) - x[i] = dx * i; - - HGPhaseFunction_2D(f, x, n, g); - - w[0][0][region] = 2.0 / 6.0 * (2.0 * f[0] + f[1]) * dx; - w[0][n - 1][region] = 2.0 / 6.0 * (f[n - 2] + 2.0 * f[n - 1]) * dx; - - for (i = 1; i < n - 1; i++) - { - w[0][i][region] = 1.0 / 6.0 * (f[i - 1] + 2.0 * f[i]) * dx + 1.0 / 6.0 * (2.0 * f[i] + f[i + 1]) * dx; - w[0][2 * n - 2 - i][region] = w[0][i][region]; - } - - norm = 0; - for (i = 0; i < nAngle; i++) - norm += w[0][i][region]; - - //Normalize - for (i = 0; i < nAngle; i++) - w[0][i][region] = w[0][i][region] / norm; - - - - for (i = 1; i < nAngle; i++) - for (j = 0; j < nAngle; j++) - { - nj = (j + i) % nAngle; - if (nj == 0) - { - nj = 0; - } - w[i][nj][region] = w[0][j][region]; - } - - for (i = 0; i < nAngle; i++) - { - theta[i][2] = 2 * 3.14159265 * i / nAngle; - theta[i][0] = Math.Cos(theta[i][2]); - theta[i][1] = Math.Sin(theta[i][2]); - } - } - - /// - /// This function describes phase function or scattering kernel and is set to be Henyey-Greenstein phase function - /// - /// function output - /// angular array - /// number of angles - /// anisotropy factor - - public static void HGPhaseFunction_2D(double[] f, double[] dAng, int n, double g) - { - int i; - - for (i = 0; i < n; i++) - f[i] = 1 / (2 * 3.14159265) * (1 - g * g) / (1 + g * g - 2 * g * Math.Cos(dAng[i])); - } - - /// - /// compute sweep ordering - /// - /// spatial mesh - /// angular mesh - /// spatial mesh level - /// angular mesh level - public static void SweepOrdering(ref SpatialMesh[] smesh, AngularMesh[] amesh, int sMeshLevel, int aMeshLevel) - { - int i, j, k; - double[][] centOfTri; - double[] temp; - int[] idxtemp; - for (i = 0; i <= sMeshLevel; i++) - { - centOfTri = new double[smesh[i].Nt][]; - for (j = 0; j < smesh[i].Nt; j++) - { - centOfTri[j] = new double[2]; - for (k = 0; k < 2; k++) - centOfTri[j][k] = ( - smesh[i].P[smesh[i].T[j][0]][k] - + smesh[i].P[smesh[i].T[j][1]][k] - + smesh[i].P[smesh[i].T[j][2]][k]) / 3; - - } - - temp = new double [smesh[i].Nt]; - idxtemp = new int[smesh[i].Nt]; - smesh[i].So = new int[amesh[aMeshLevel].Ns][]; - for (j = 0; j < amesh[aMeshLevel].Ns; j++) - smesh[i].So[j] = new int[smesh[i].Nt]; - - for (j = 0; j < amesh[aMeshLevel].Ns; j++) - { - for (k = 0; k < smesh[i].Nt; k++) - { - temp[k] = Math.Cos(amesh[aMeshLevel].Ang[j][2]) * centOfTri[k][0] + Math.Sin(amesh[aMeshLevel].Ang[j][2]) * centOfTri[k][1]; - idxtemp[k] = k; - } - QuickSort(temp, idxtemp, 0, smesh[i].Nt - 1); - for (k = 0; k < smesh[i].Nt; k++) - smesh[i].So[j][k] = idxtemp[k]; - } - } - } - - /// - /// convert triangular mesh to a square mesh - /// - /// spatial mesh - /// x coordinates - /// y coordinates - /// grid location - /// average intensity of each traingle - /// total points - public static void SquareTriMeshToGrid(ref SpatialMesh smesh, ref double[] x, ref double[] y, ref double[][] uxy, double[] inten, int nxy) - { - int i, j, k; - int np, nt; - double dx, dy; - int[][] tn; - double[][] a12; - double[][] a13; - double[] a2 = new double[2]; - double[] a3 = new double[2]; - double[] b2 = new double[2]; - double[] b3 = new double[2]; - double d2, d3; - double[] r1p = new double[2]; - - double xmin, xmax, zmin, zmax; - double temp; - double tiny = 2.2204e-12; - - double[] xminArray; - double[] xmaxArray; - double[] zminArray; - double[] zmaxArray; - - int[] idxxminArray; - int[] idxxmaxArray; - int[] idxzminArray; - int[] idxzmaxArray; - - // Assign large number for min and small number for max - xmin = 1e10; xmax = -1e10; zmin = 1e10; zmax = -1e10; - - np = smesh.Np; - nt = smesh.Nt; - - for (i = 0; i < np; i++) - { - xmin = Math.Min(smesh.P[i][0], xmin); - xmax = Math.Max(smesh.P[i][0], xmax); - zmin = Math.Min(smesh.P[i][1], zmin); - zmax = Math.Max(smesh.P[i][1], zmax); - } - - dx = (xmax - xmin) / (nxy - 1); - dy = (zmax - zmin) / (nxy - 1); - - for (i = 0; i < nxy; i++) - { - x[i] = xmin + i * dx; - y[i] = zmin + i * dy; - } - - xminArray = new double[nt]; - xmaxArray = new double[nt]; - zminArray = new double[nt]; - zmaxArray = new double[nt]; - - idxxminArray = new int[nt]; - idxxmaxArray = new int[nt]; - idxzminArray = new int[nt]; - idxzmaxArray = new int[nt]; - - for (i = 0; i < nt; i++) - { - xmin = 1e10; xmax = -1e10; zmin = 1e10; zmax = -1e10; - - xmin = Math.Min(smesh.P[smesh.T[i][0]][0], xmin); - xmin = Math.Min(smesh.P[smesh.T[i][1]][0], xmin); - xmin = Math.Min(smesh.P[smesh.T[i][2]][0], xmin); - - xmax = Math.Max(smesh.P[smesh.T[i][0]][0], xmax); - xmax = Math.Max(smesh.P[smesh.T[i][1]][0], xmax); - xmax = Math.Max(smesh.P[smesh.T[i][2]][0], xmax); - - zmin = Math.Min(smesh.P[smesh.T[i][0]][1], zmin); - zmin = Math.Min(smesh.P[smesh.T[i][1]][1], zmin); - zmin = Math.Min(smesh.P[smesh.T[i][2]][1], zmin); - - zmax = Math.Max(smesh.P[smesh.T[i][0]][1], zmax); - zmax = Math.Max(smesh.P[smesh.T[i][1]][1], zmax); - zmax = Math.Max(smesh.P[smesh.T[i][2]][1], zmax); - - xminArray[i] = xmin; - xmaxArray[i] = xmax; - zminArray[i] = zmin; - zmaxArray[i] = zmax; - - idxxminArray[i] = i; - idxxmaxArray[i] = i; - idxzminArray[i] = i; - idxzmaxArray[i] = i; - } - - QuickSort(xminArray, idxxminArray, 0, nt - 1); - QuickSort(xmaxArray, idxxmaxArray, 0, nt - 1); - QuickSort(zminArray, idxzminArray, 0, nt - 1); - QuickSort(zmaxArray, idxzmaxArray, 0, nt - 1); - - j = 0; - for (i = 0; i < nt; i++) - { - if (j < nxy) - { - while (x[j] < xminArray[i]) - { - j++; - if (j >= nxy) - break; - } - } - xminArray[i] = j; - } - - j = nxy - 1; - for (i = nt - 1; i >= 0; i--) - { - if (j >= 0) - { - while (x[j] > xmaxArray[i]) - { - j--; - if (j < 0) - break; - } - } - xmaxArray[i] = j; - } - - - j = 0; - for (i = 0; i < nt; i++) - { - if (j < nxy) - { - while (y[j] < zminArray[i]) - { - j++; - if (j >= nxy) - break; - } - } - zminArray[i] = j; - } - - j = nxy - 1; - for (i = nt - 1; i >= 0; i--) - { - if (j >= 0) - { - while (y[j] > zmaxArray[i]) - { - j--; - if (j < 0) - break; - } - } - zmaxArray[i] = j; - } - - RearrangeArray(ref xminArray, idxxminArray, nt); - RearrangeArray(ref xmaxArray, idxxmaxArray, nt); - RearrangeArray(ref zminArray, idxzminArray, nt); - RearrangeArray(ref zmaxArray, idxzmaxArray, nt); - - - tn = new int[nxy][]; - a12 = new double[nxy][]; - a13 = new double[nxy][]; - for (i = 0; i < nxy; i++) - { - tn[i] = new int[nxy]; - a12[i] = new double[nxy]; - a13[i] = new double[nxy]; - } - - //Set tn to a non positive number - for (i = 0; i < nxy; i++) - for (j = 0; j < nxy; j++) - { - tn[i][j] = -1; - } - - - for (i = 0; i < nt; i++) - { - if ((xminArray[i] <= xmaxArray[i]) && (zminArray[i] <= zmaxArray[i])) - { - for (j = 0; j < 2; j++) - { - a2[j] = smesh.P[smesh.T[i][1]][j] - smesh.P[smesh.T[i][0]][j]; - a3[j] = smesh.P[smesh.T[i][2]][j] - smesh.P[smesh.T[i][0]][j]; - } - - temp = a2[0] * a3[1] - a2[1] * a3[0]; - b2[0] = a3[1] / temp; - b2[1] = -a3[0] / temp; - b3[0] = -a2[1] / temp; - b3[1] = a2[0] / temp; - - for (j = (int)xminArray[i]; j <= (int)xmaxArray[i]; j++) - { - for (k = (int)zminArray[i]; k <= (int)zmaxArray[i]; k++) - { - if (tn[k][j] == -1) - { - r1p[0] = x[j] - smesh.P[smesh.T[i][0]][0]; - r1p[1] = y[k] - smesh.P[smesh.T[i][0]][1]; - d2 = b2[0] * r1p[0] + b2[1] * r1p[1]; - if ((d2 >= -tiny) && (d2 <= 1 + tiny)) - { - d3 = b3[0] * r1p[0] + b3[1] * r1p[1]; - if ((d3 >= -tiny) && (d2 + d3 <= 1 + tiny)) - { - tn[k][j] = i; - a12[k][j] = d2; - a13[k][j] = d3; - } - } - } - } - } - } - } - - - double tt1, tt2, tt3; - - for (i = 0; i < nxy; i++) - for (j = 0; j < nxy; j++) - { - tt1 = (1 - a12[i][j] - a13[i][j]) * inten[(int)smesh.T[tn[i][j]][0]]; - tt2 = a12[i][j] * inten[(int)smesh.T[tn[i][j]][1]]; - tt3 = a13[i][j] * inten[(int)smesh.T[tn[i][j]][2]]; - uxy[i][j] = tt1 + tt2 + tt3; - - } - } - - private static void QuickSort(double[] a, int[] idx, int n_left, int n_right) - { - int i, j; - double pivot; - - if (n_left >= n_right) - return; - i = n_left; - j = n_right; - - // Choose the last element as pivot - pivot = a[j]; - while (i < j) - { - while (i < j && a[i] <= pivot) i++; - while (i < j && a[j] >= pivot) j--; - if (i < j) - { - Swap(ref a[i], ref a[j]); - Swap(ref idx[i], ref idx[j]); - } - } - if (i != n_right) - { - Swap(ref a[i], ref a[n_right]); - Swap(ref idx[i], ref idx[n_right]); - } - - // sort sub-array recursively - QuickSort(a, idx, n_left, i - 1); - QuickSort(a, idx, i + 1, n_right); - } - - private static void Swap(ref double x, ref double y) - { - double temp = y; - y = x; - x = temp; - } - - private static void Swap(ref int x, ref int y) - { - int temp = y; - y = x; - x = temp; - } - - /// - /// Rearrange an array based on its index array - /// - /// data array - /// index of data array - /// number of elements in the array - private static void RearrangeArray(ref double[] a, int[] idx, int n) - { - int i; - double[] temp = new double[n]; - - for (i = 0; i < n; i++) - temp[idx[i]] = a[i]; - a = temp; - } - - - /// - /// Create an angular mesh - /// - /// - /// - /// - public static void CreateAnglularMesh(ref AngularMesh[] amesh, int aLevel, MultiEllipsoidTissueInput tissueInput) - { - int i, j, k, l; - int nLayerRegions = tissueInput.LayerRegions.Length; - int nInclusionRegions = tissueInput.EllipsoidRegions.Length; - int nRegions = nLayerRegions + nInclusionRegions - 2; - - - for (i = 0; i <= aLevel; i++) - { - amesh[i].Ns = (int) Math.Pow(2.0, (double) (i + 1)); - amesh[i].Ang = new double[amesh[i].Ns][]; - amesh[i].W = new double[amesh[i].Ns][][]; - for (j = 0; j < amesh[i].Ns; j++) - { - amesh[i].Ang[j] = new double[3]; - amesh[i].W[j] = new double[amesh[i].Ns][]; - for (k = 0; k < amesh[i].Ns; k++) - amesh[i].W[j][k] = new double[nRegions]; - } - - for (j = 0; j < amesh[i].Ns; j++) - { - for (l = 0; l < nLayerRegions - 2; l++) - Weight_2D(amesh[i].W, amesh[i].Ang, amesh[i].Ns, tissueInput.LayerRegions[l + 1].RegionOP.G, l); - - for (l = 0; l < nInclusionRegions; l++) - Weight_2D(amesh[i].W, amesh[i].Ang, amesh[i].Ns, tissueInput.EllipsoidRegions[l].RegionOP.G, - l + nLayerRegions - 2); - } - } - } - - /// - /// Assign region values to smesh[].Region variable - /// - /// spatial mesh - /// number of nodes - /// tissue input - public static void AssignRegions(ref SpatialMesh[] smesh, int np, MultiEllipsoidTissueInput tissueInput) - { - int i, j, k; - double x, z; - - int nLayers = tissueInput.LayerRegions.Length; - int nInclusions = tissueInput.EllipsoidRegions.Length; - - for (i = 0; i <= np; i++) - for (j = 0; j < smesh[i].Np; j++) - { - x = smesh[i].P[j][0]; - z = smesh[i].P[j][1]; - - for (k = 1; k < nLayers - 1; k++) - { - if ((z >= ((LayerTissueRegion)tissueInput.LayerRegions[k]).ZRange.Start) && - (z <= ((LayerTissueRegion)tissueInput.LayerRegions[k]).ZRange.Stop)) - smesh[i].Region[j] = k - 1; - - } - - for (k = 0; k < nInclusions; k++) - { - double dx = ((EllipsoidTissueRegion)tissueInput.EllipsoidRegions[k]).Dx; - double dz = ((EllipsoidTissueRegion)tissueInput.EllipsoidRegions[k]).Dz; - - Position center = ((EllipsoidTissueRegion)tissueInput.EllipsoidRegions[k]).Center; - - double func = ((x - center.X) * (x - center.X) / (dx * dx)) + ((z - center.Z) * (z - center.Z) / (dz * dz)); - if (func <= 1) - smesh[i].Region[j] = nLayers - 2 + k; - - } - } - } - - /// - /// Set Mus in nodes - /// - /// mua - /// spatial mesh - /// simulation input - public static void SetMua(ref double[][][] ua, SpatialMesh[] smesh, SimulationInput input) - { - int j, k, l; - int nt = smesh[input.MeshDataInput.SMeshLevel].Nt; - var tissueInput = ((MultiEllipsoidTissueInput) input.TissueInput); - int nLayerRegions = tissueInput.LayerRegions.Length; - int nInclusionRegions = tissueInput.EllipsoidRegions.Length; - int sMeshLevel = input.MeshDataInput.SMeshLevel; - - ua[sMeshLevel] = new double[nt][]; - for (j = 0; j < nt; j++) - { - ua[sMeshLevel][j] = new double[3]; - for (k = 0; k < 3; k++) - { - for (l = 0; l < nLayerRegions - 2; l++) - { - if (smesh[sMeshLevel].Region[smesh[sMeshLevel].T[j][0]] == l) - ua[sMeshLevel][j][k] = tissueInput.LayerRegions[l+1].RegionOP.Mua; - } - - for (l = 0; l < nInclusionRegions; l++) - { - if (smesh[sMeshLevel].Region[smesh[sMeshLevel].T[j][0]] == l + nLayerRegions - 2) - ua[sMeshLevel][j][k] = tissueInput.EllipsoidRegions[l].RegionOP.Mua; - } - } - } - } - - /// - /// Set Mus in nodes - /// - /// mus - /// spatial mesh - /// simulation input - public static void SetMus(ref double[][][] us, SpatialMesh[] smesh, SimulationInput input) - { - int j, k, l; - int nt = smesh[input.MeshDataInput.SMeshLevel].Nt; - var tissueInput = ((MultiEllipsoidTissueInput)input.TissueInput); - int nLayerRegions = tissueInput.LayerRegions.Length; - int nInclusionRegions = tissueInput.EllipsoidRegions.Length; - int sMeshLevel = input.MeshDataInput.SMeshLevel; - - us[sMeshLevel] = new double[nt][]; - for (j = 0; j < nt; j++) - { - us[sMeshLevel][j] = new double[3]; - for (k = 0; k < 3; k++) - { - for (l = 0; l < nLayerRegions - 2; l++) - { - if (smesh[sMeshLevel].Region[smesh[sMeshLevel].T[j][0]] == l) - us[sMeshLevel][j][k] = tissueInput.LayerRegions[l + 1].RegionOP.Mus; - } - - for (l = 0; l < nInclusionRegions; l++) - { - if (smesh[sMeshLevel].Region[smesh[sMeshLevel].T[j][0]] == l + nLayerRegions - 2) - us[sMeshLevel][j][k] = tissueInput.EllipsoidRegions[l].RegionOP.Mus; - } - } - } - } - - - /// - /// Create a squarte mesh for given spatial mesh level - /// - /// spatial mesh - /// spatial mesh levels - /// length - public static void CreateSquareMesh(ref SpatialMesh[] smesh, int sMeshLevel, double length) - { - int i; - int np, nt, ne; - - //SQUARE MESH - np = 5; - ne = 4; - nt = 4; - - double[][] pts = new double[np][]; - for (i = 0; i < np; i++) - pts[i] = new double[4]; - - int[][] edge = new int[ne][]; - for (i = 0; i < ne; i++) - edge[i] = new int[4]; - - int[][] tri = new int[nt][]; - for (i = 0; i < nt; i++) - tri[i] = new int[3]; - - //create the basic square mesh - BasicSquareMesh(pts, edge, tri, length); - - AssignSpatialMesh(ref smesh, pts, edge, tri, np, ne, nt, 0); - - //string str = "PET" + 0; - //str = str + ".txt"; - //WritePTEData(str, pts, edge, tri, np, ne, nt); - - if (sMeshLevel > 0) - CreateMultigrid(ref smesh, pts, edge, tri, np, ne, nt, sMeshLevel); - } - - /// - /// Define Basic Square Mesh - /// - /// points data - /// edge data - /// triangle data - /// length - private static void BasicSquareMesh( - double[][] pts, - int[][] edge, - int[][] tri, - double length) - { - //Boundaries: For X: x =-0.5; x =0.5 & For Z: z=0; z=1; - pts[0][0] = -0.5 *length; pts[0][1] = 0; pts[0][2] = 0; pts[0][3] = 1; - pts[1][0] = 0.5 * length; pts[1][1] = 0; pts[1][2] = 1; pts[1][3] = 1; - pts[2][0] = 0.5 * length; pts[2][1] = length; pts[2][2] = 2; pts[2][3] = 1; - pts[3][0] = -0.5 *length; pts[3][1] = length; pts[3][2] = 3; pts[3][3] = 1; - pts[4][0] = 0; pts[4][1] = 0.5 * length; pts[4][2] = 4; pts[4][3] = 0; - - edge[0][0] = -1; edge[0][1] = (int)pts[0][2]; edge[0][2] = (int)pts[1][2]; edge[0][3] = -1; - edge[1][0] = -1; edge[1][1] = (int)pts[1][2]; edge[1][2] = (int)pts[2][2]; edge[1][3] = -1; - edge[2][0] = -1; edge[2][1] = (int)pts[2][2]; edge[2][2] = (int)pts[3][2]; edge[2][3] = -1; - edge[3][0] = -1; edge[3][1] = (int)pts[3][2]; edge[3][2] = (int)pts[0][2]; edge[3][3] = -1; - - tri[0][0] = (int)pts[0][2]; tri[0][1] = (int)pts[1][2]; tri[0][2] = (int)pts[4][2]; - tri[1][0] = (int)pts[1][2]; tri[1][1] = (int)pts[2][2]; tri[1][2] = (int)pts[4][2]; - tri[2][0] = (int)pts[2][2]; tri[2][1] = (int)pts[3][2]; tri[2][2] = (int)pts[4][2]; - tri[3][0] = (int)pts[3][2]; tri[3][1] = (int)pts[0][2]; tri[3][2] = (int)pts[4][2]; - } - - /// - /// Create Multigrid based on spatial mesh level - /// - /// - /// points data - /// edge data - /// triangle data - /// number of points - /// number of boundary edges - /// number of triangles - /// spatial mesh levels - public static void CreateMultigrid(ref SpatialMesh[] smesh, double[][] p, int[][] e, int[][] t, int np, int ne, int nt, int sMeshLevel) - { - int i, j, nt2, ne2, np2, npt, trint; - - //Input arrays - double[][] oldp = p; - int[][] oldt = t; - int[][] olde = e; - - //char[] str[80]; - //char[] level[4]; - - for (j = 1; j <= sMeshLevel; j++) - { - //npt gets the total number of point after sub division - npt = 0; - - trint = 3 * nt; - - double[][] ptemp = new double[trint][]; - for (i = 0; i < trint; i++) - ptemp[i] = new double[6]; - - FindPTemp(ptemp, oldp, oldt, nt); - ReindexingPTemp(ptemp, ref npt, trint, np); - - np2 = npt; - ne2 = 2 * ne; - nt2 = 4 * nt; - - double[][] newp = new double[np2][]; - for (i = 0; i < np2; i++) - newp[i] = new double[4]; - - int[][] newe = new int[ne2][]; - for (i = 0; i < ne2; i++) - newe[i] = new int[4]; - - int[][] newt = new int[nt2][]; - for (i = 0; i < nt2; i++) - newt[i] = new int[3]; - - NewPET(newp, oldp, newe, newt, oldt, ptemp, np, nt, trint); - - //Update np, nt and ne - np = np2; - ne = ne2; - nt = nt2; - - //Assign output arrays to input arrays - oldp = newp; - olde = newe; - oldt = newt; - - AssignSpatialMesh(ref smesh, oldp, olde, oldt, np, ne, nt, j); - } - - for (j = 0; j <= sMeshLevel; j++) - smesh[j].Region = new int[smesh[j].Np]; - - - - } - - private static void FindPTemp(double[][] ptemp, double[][] p, int[][] t, int nt) - { - int i; - int p0, p1, p2; - - for (i = 0; i < nt; i++) - { - //x cordinates - ptemp[3 * i][0] = 0.5 * (p[t[i][0]][0] + p[t[i][1]][0]); - ptemp[3 * i + 1][0] = 0.5 * (p[t[i][0]][0] + p[t[i][2]][0]); - ptemp[3 * i + 2][0] = 0.5 * (p[t[i][1]][0] + p[t[i][2]][0]); - - //y cordinates - ptemp[3 * i][1] = 0.5 * (p[t[i][0]][1] + p[t[i][1]][1]); - ptemp[3 * i + 1][1] = 0.5 * (p[t[i][0]][1] + p[t[i][2]][1]); - ptemp[3 * i + 2][1] = 0.5 * (p[t[i][1]][1] + p[t[i][2]][1]); - - //P index - initialize to negative values - ptemp[3 * i][2] = -1; - ptemp[3 * i + 1][2] = -1; - ptemp[3 * i + 2][2] = -1; - - //find the edge vertex of a triangle - p0 = (int)p[t[i][0]][3]; - p1 = (int)p[t[i][1]][3]; - p2 = (int)p[t[i][2]][3]; - - if (p0 + p1 == 2) - { - ptemp[3 * i][3] = 1; - ptemp[3 * i][4] = p[t[i][0]][2]; - ptemp[3 * i][5] = p[t[i][1]][2]; - } - else if (p0 + p2 == 2) - { - ptemp[3 * i + 1][3] = 1; - ptemp[3 * i + 1][4] = p[t[i][0]][2]; - ptemp[3 * i + 1][5] = p[t[i][2]][2]; - } - else if (p1 + p2 == 2) - { - ptemp[3 * i + 2][3] = 1; - ptemp[3 * i + 2][4] = p[t[i][1]][2]; - ptemp[3 * i + 2][5] = p[t[i][2]][2]; - } - } - - } - - private static void ReindexingPTemp(double[][] ptemp, ref int count, int trint, int np) - { - int i, j; - - count = np; - for (i = 0; i < trint - 1; i++) - { - for (j = i + 1; j < trint; j++) - { - if (ptemp[i][0] == ptemp[j][0]) - if (ptemp[i][1] == ptemp[j][1]) - { - ptemp[j][2] = count; - ptemp[i][2] = count; - count++; - } - } - } - for (i = 0; i < trint; i++) - { - if (ptemp[i][2] < 0) - { - ptemp[i][2] = count; - count++; - } - } - } - - /// - /// Assign P E T - /// - private static void NewPET(double[][] newp, double[][] oldp, int[][] newe, int[][] newt, int[][] oldt, double[][] ptemp, int np, int nt, int trint) - { - int i, count; - - //Assign P - for (i = 0; i < np; i++) - { - count = (int)oldp[i][2]; - newp[count][0] = oldp[i][0]; - newp[count][1] = oldp[i][1]; - newp[count][2] = oldp[i][2]; - newp[count][3] = oldp[i][3]; - } - - for (i = 0; i < trint; i++) - { - count = (int)ptemp[i][2]; - newp[count][0] = ptemp[i][0]; - newp[count][1] = ptemp[i][1]; - newp[count][2] = ptemp[i][2]; - newp[count][3] = ptemp[i][3]; - } - - - //Assign E - count = 0; - for (i = 0; i < trint; i++) - { - if (ptemp[i][3] == 1) - { - newe[count][0] = -1; - newe[count][1] = (int)ptemp[i][4]; - newe[count][2] = (int)ptemp[i][2]; - newe[count][3] = -1; - count++; - newe[count][0] = -1; - newe[count][1] = (int)ptemp[i][2]; - newe[count][2] = (int)ptemp[i][5]; - newe[count][3] = -1; - count++; - } - } - - - - double x0, x1, x2, x3, x4, x5; - double y0, y1, y2, y3, y4, y5; - - double temp; - - //Assign T - for (i = 0; i < nt; i++) - { - - x0 = oldp[oldt[i][0]][0]; - y0 = oldp[oldt[i][0]][1]; - - x1 = oldp[oldt[i][1]][0]; - y1 = oldp[oldt[i][1]][1]; - - x2 = oldp[oldt[i][2]][0]; - y2 = oldp[oldt[i][2]][1]; - - x3 = ptemp[3 * i][0]; - y3 = ptemp[3 * i][1]; - - x4 = ptemp[3 * i + 1][0]; - y4 = ptemp[3 * i + 1][1]; - - x5 = ptemp[3 * i + 2][0]; - y5 = ptemp[3 * i + 2][1]; - - //Nodes in the triangle should be arranged counter clock wise - // http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm - //(x1-x0)*(y2-y0)-(x2-x0)*(y1-y0) > 0 counter clockwise - - //Triangle 1 - temp = (x3 - x0) * (y4 - y0) - (x4 - x0) * (y3 - y0); - if (temp > 0) - { - newt[4 * i][0] = oldt[i][0]; - newt[4 * i][1] = (int)ptemp[3 * i][2]; - newt[4 * i][2] = (int)ptemp[3 * i + 1][2]; - } - else - { - newt[4 * i][0] = oldt[i][0]; - newt[4 * i][1] = (int)ptemp[3 * i + 1][2]; - newt[4 * i][2] = (int)ptemp[3 * i][2]; - } - - //Triangle 2 - temp = (x3 - x1) * (y5 - y1) - (x5 - x1) * (y3 - y1); - if (temp > 0) - { - newt[4 * i + 1][0] = oldt[i][1]; - newt[4 * i + 1][1] = (int)ptemp[3 * i][2]; - newt[4 * i + 1][2] = (int)ptemp[3 * i + 2][2]; - } - else - { - newt[4 * i + 1][0] = oldt[i][1]; - newt[4 * i + 1][1] = (int)ptemp[3 * i + 2][2]; - newt[4 * i + 1][2] = (int)ptemp[3 * i][2]; - } - - //Triangle 3 - temp = (x4 - x2) * (y5 - y2) - (x5 - x2) * (y4 - y2); - if (temp > 0) - { - newt[4 * i + 2][0] = oldt[i][2]; - newt[4 * i + 2][1] = (int)ptemp[3 * i + 1][2]; - newt[4 * i + 2][2] = (int)ptemp[3 * i + 2][2]; - } - else - { - newt[4 * i + 2][0] = oldt[i][2]; - newt[4 * i + 2][1] = (int)ptemp[3 * i + 2][2]; - newt[4 * i + 2][2] = (int)ptemp[3 * i + 1][2]; - } - - //Triangle 4 (Inner Triangle) - temp = (x4 - x3) * (y5 - y3) - (x5 - x3) * (y4 - y3); - if (temp > 0) - { - newt[4 * i + 3][0] = (int)ptemp[3 * i + 1][2]; - newt[4 * i + 3][1] = (int)ptemp[3 * i + 2][2]; - newt[4 * i + 3][2] = (int)ptemp[3 * i][2]; - } - else - { - newt[4 * i + 3][0] = (int)ptemp[3 * i + 1][2]; - newt[4 * i + 3][1] = (int)ptemp[3 * i][2]; - newt[4 * i + 3][2] = (int)ptemp[3 * i + 2][2]; - } - } - } - - - private static void AssignSpatialMesh(ref SpatialMesh[] smesh, double[][] p, int[][] e, int[][] t, int np, int ne, int nt, int level) - { - int i, j; - - smesh[level].Np = np; - smesh[level].Ne = ne; - smesh[level].Nt = nt; - - smesh[level].P = new double[np][]; - for (i = 0; i < np; i++) - { - smesh[level].P[i] = new double[2]; - for (j = 0; j < 2; j++) - smesh[level].P[i][j] = p[i][j]; - } - - smesh[level].T = new int[nt][]; - for (i = 0; i < nt; i++) - { - smesh[level].T[i] = new int[3]; - for (j = 0; j < 3; j++) - smesh[level].T[i][j] = t[i][j]; - } - - - smesh[level].E = new int[ne][]; - for (i = 0; i < ne; i++) - { - smesh[level].E[i] = new int[4]; - for (j = 0; j < 4; j++) - smesh[level].E[i][j] = e[i][j]; - } - - } - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/MultiGridCycle.cs b/src/Vts.FemModeling/MGRTE/2D/MultiGridCycle.cs deleted file mode 100644 index 81fa327ed..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/MultiGridCycle.cs +++ /dev/null @@ -1,1668 +0,0 @@ -using System; -using Vts.FemModeling.MGRTE._2D.DataStructures; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// This class contains numerical algorithms required to solve MG-RTE. - /// - public class MultiGridCycle - { - - private double Pi = GlobalConstants.Pi; - - public void Boundary(int ne, int nt, int[][] t, int[][] p2, double[][] p, int[][] e, int[][] e2, int[][] so2, double[][] n, int[] ori) - // Purpose: this function is to compute "e", "e2", "so2", "n", "sn" and "ori" with the data structure as follow: - // - // e[ne][4]: for the 2nd index,1st is the corresponding triangle containing boundary, - // 2nd and 3rd are global nodal index of boundary nodes, and 4th - // is global index of the other non-boundary node in the triangle. - // - // e2[ne][2]: the local order of the boundary nodes in the corresponding triangle with - // 1st entry for node 1 (e[ne][1]) and 2nd entry for node 2 (e[ne][2]). - // - // so2[nt][3]: the corresponding triangle and edge of the boundary - // Example: if boundary "10" is edge "2" of triangle "20", - // then so2[20][2]=10, and "-1" otherwise. - // - // n[ne][2]: x and y coordinate of the outgoing normal. - // - // ori[ne]: orientation of the boundary edge w.r.t. the triangle, which is useful at ONLY one place in "bc_reflection". - { - int i, j, k, tri; - double nx, ny, temp; - - for (i = 0; i < nt; i++) - { - for (j = 0; j < 3; j++) - { - so2[i][j] = -1; - } - } - for (i = 0; i < ne; i++) - { - tri = FindTri2(p2[e[i][1]], p2[e[i][2]]); - e[i][0] = tri; - - // find the local index for boundary nodes in the corresponding triangle - for (j = 0; j < 2; j++) - { - for (k = 0; k < 3; k++) - { - if (e[i][j + 1] == t[tri][k]) - { - e2[i][j] = k; goto stop; - } - } - stop: ; - } - // find the non-boundary node and the boundary index in the corresponding triangle - for (j = 0; j < 3; j++) - { - if (t[tri][j] != e[i][1] && t[tri][j] != e[i][2]) - { - e[i][3] = t[tri][j]; //non-edge node - so2[tri][j] = i; //boundary index - goto stop2; - } - } - stop2: ; - - // The following is to compute outgoing vector "n" on the boundary - nx = -(p[e[i][1]][1] - p[e[i][2]][1]); - ny = p[e[i][1]][0] - p[e[i][2]][0]; - ori[i] = 1; - if (nx * (p[e[i][3]][0] - p[e[i][2]][0]) + ny * (p[e[i][3]][1] - p[e[i][2]][1]) > 0)// the normal is incoming. - { - nx = -nx; ny = -ny; ori[i] = 0; - } - temp = Math.Sqrt(nx * nx + ny * ny); - n[i][0] = nx / temp; n[i][1] = ny / temp; - } - } - - - public void BoundReflection(int ns, double[][] theta, SpatialMesh smesh, double nTissue, double nExt, BoundaryCoupling b) - // Purpose: this fucntion is to find the coupling relation between directions on the boundary - // due to reflection and refraction in the presence of refraction index mismatch at the boundary. - // For the data structure of "b", see "struct boundarycoupling" in "solver". - { - int i, j; - int[] sign = new int[2] { -1, 1 }; - double dx, dy, sn, dtheta = 2 * Pi / ns, ratio_reflection, ratio_refraction; - double[] temp = new double[2]; - double theta0, theta_i, theta_m, theta_m2; - - for (i = 0; i < smesh.Ne; i++) - { - dx = smesh.P[smesh.E[i][1]][0] - smesh.P[smesh.E[i][2]][0]; - dy = smesh.P[smesh.E[i][1]][1] - smesh.P[smesh.E[i][2]][1]; - if (smesh.Ori[i] == 0) - { - dx = -dx; dy = -dy; // to make sure that (dx,dy) goes clockwisely. - } - for (j = 0; j < ns; j++) - { - sn = theta[j][0] * smesh.N[i][0] + theta[j][1] * smesh.N[i][1]; // "s" dot "n" for the angle "s" - if (sn < 0) - { - theta0 = Pi - Math.Acos(sn); - ratio_reflection = Reflection(theta0, nTissue, nExt); - Refraction(temp, theta0, nExt, nTissue); - ratio_refraction = temp[0]; theta_i = temp[1]; - - if (theta[j][0] * dx + theta[j][1] * dy > 0) // the ONLY place for clockwise (dx,dy) - { - theta_m = Mod2pi(theta[j][2] + (Pi - 2 * theta0)); - theta_m2 = Mod2pi(theta[j][2] + (theta_i - theta0)); - } - else - { - theta_m = Mod2pi(theta[j][2] - (Pi - 2 * theta0)); - theta_m2 = Mod2pi(theta[j][2] - (theta_i - theta0)); - } - // contribution to incoming flux through internal reflection - MathFunctions.Intepolation_a(theta_m, dtheta, ns, b.Ri[i][j], b.Ri2[i][j], ratio_reflection); - // contribution from boundary source to incoming flux through refraction - MathFunctions.Intepolation_a(theta_m2, dtheta, ns, b.Si[i][j], b.Si2[i][j], ratio_refraction); - } - else - { - theta0 = Math.Acos(sn); - ratio_reflection = Reflection(theta0, nExt, nTissue); - Refraction(temp, theta0, nTissue, nExt); - ratio_refraction = temp[0]; theta_i = temp[1]; - - if (theta[j][0] * dx + theta[j][1] * dy > 0) // the ONLY place for clockwise (dx,dy) - { - theta_m = Mod2pi(theta[j][2] - (Pi - 2 * theta0)); - theta_m2 = Mod2pi(theta[j][2] - (theta_i - theta0)); - } - else - { - theta_m = Mod2pi(theta[j][2] + (Pi - 2 * theta0)); - theta_m2 = Mod2pi(theta[j][2] + (theta_i - theta0)); - } - // contribution from boundary source to outgoing flux through reflection - MathFunctions.Intepolation_a(theta_m, dtheta, ns, b.Ro[i][j], b.Ro2[i][j], ratio_reflection); - // contribution to outgoing flux after refraction - MathFunctions.Intepolation_a(theta_m2, dtheta, ns, b.So[i][j], b.So2[i][j], ratio_refraction); - } - } - } - } - - - public void CtoF(int nt_c, int ns_c, double[][][] flux, double[][][] dcflux, int[][] smap, double[][][][] cf) - // Purpose: this function describes coarse-to-fine operator simultaneouly in angle and space, - // i.e., the combination of "ctof_a" and "ctof_s". - { - int i, j, k, m, n; - for (i = 0; i < ns_c - 1; i++) - { - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - flux[2 * i][smap[j][m]][n] += dcflux[i][j][k] * cf[j][k][m - 1][n]; - flux[2 * i + 1][smap[j][m]][n] += 0.5 * (dcflux[i][j][k] + dcflux[i + 1][j][k]) * cf[j][k][m - 1][n]; - } - } - } - } - } - i = ns_c - 1; - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - flux[2 * i][smap[j][m]][n] += dcflux[i][j][k] * cf[j][k][m - 1][n]; - flux[2 * i + 1][smap[j][m]][n] += 0.5 * (dcflux[i][j][k] + dcflux[0][j][k]) * cf[j][k][m - 1][n]; - } - } - } - } - } - - - public void CtoF_a(int nt, int ns_c, double[][][] flux, double[][][] dcflux) - // Purpose: this function describes angular coarse-to-fine operator by linear interpolation. - { - int i, j, k; - for (i = 0; i < nt; i++) - { - for (j = 0; j < ns_c - 1; j++) - { - for (k = 0; k < 3; k++) - { - flux[2 * j][i][k] += dcflux[j][i][k]; - flux[2 * j + 1][i][k] += (dcflux[j][i][k] + dcflux[j + 1][i][k]) / 2; - } - } - j = ns_c - 1; - for (k = 0; k < 3; k++) - { - flux[2 * j][i][k] += dcflux[j][i][k]; - flux[2 * j + 1][i][k] += (dcflux[j][i][k] + dcflux[0][i][k]) / 2; - } - } - } - - - public void CtoF_s(int nt_c, int ns, double[][][] flux, double[][][] dcflux, int[][] smap, double[][][][] cf) - // Purpose: this function describes spatial coarse-to-fine operator by linear interpolation for piecewise linear DG. - { - int i, j, k, m, n; - for (i = 0; i < ns; i++) - { - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - flux[i][smap[j][m]][n] += dcflux[i][j][k] * cf[j][k][m - 1][n]; - } - } - } - } - } - } - - - private void CtoF_s2(int nt_c, double[][] f, double[][] cf2, int[][] smap, double[][][][] cf) - // Purpose: this function describes spatial coarse-to-fine operator by linear interpolation for piecewise linear DG. - { - int j, k, m, n; - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - f[smap[j][m]][n] = cf2[j][k] * cf[j][k][m - 1][n]; - } - } - } - } - } - - - public void Defect(AngularMesh amesh, SpatialMesh smesh, int Ns, double[][][] RHS, double[][] ua, double[][] us, double[][][] flux, - BoundaryCoupling bb, double[][][] q, double[][][] res, int vacuum) - - // Purpose: this function is to compute the residual with vacuum or reflection boundary condition. - // see "relaxation" for more details. - { - int i, j, k, m, ii, jj, bi, aMeshLevel, edge = 0; - - double[,] left = new double[3, 3]; - double[] right = new double[3]; - double[] temp = new double[3]; - double tempd; - double dettri, cosi, sini, a, b, c, source_corr; - double[,] bv = new double[3, 3] { { 1.0 / 3, 1.0 / 6, 1.0 / 6 }, { 1.0 / 6, 1.0 / 3, 1.0 / 6 }, { 1.0 / 6, 1.0 / 6, 1.0 / 3 } }; - double[,] matrix1 = new double[3, 3]; - double[,] matrix2 = new double[3, 3]; - int[,] index = new int[3, 2]; - - - index[0, 0] = 1; index[0, 1] = 2; - index[1, 0] = 2; index[1, 1] = 0; - index[2, 0] = 0; index[2, 1] = 1; - - aMeshLevel = Ns / amesh.Ns; - - if (vacuum == 0) - { - for (i = 0; i < amesh.Ns; i++) - { - bi = i * aMeshLevel; - for (j = 0; j < smesh.Nt; j++) - { - dettri = 2 * smesh.A[j]; - cosi = amesh.Ang[i][0]; sini = amesh.Ang[i][1]; - a = cosi * (smesh.P[smesh.T[j][2]][1] - smesh.P[smesh.T[j][0]][1]) + sini * (smesh.P[smesh.T[j][0]][0] - smesh.P[smesh.T[j][2]][0]); - b = cosi * (smesh.P[smesh.T[j][0]][1] - smesh.P[smesh.T[j][1]][1]) + sini * (smesh.P[smesh.T[j][1]][0] - smesh.P[smesh.T[j][0]][0]); - MatrixConvec(a, b, matrix1); - - a = ua[j][0] + us[j][0]; - b = ua[j][1] + us[j][1]; - c = ua[j][2] + us[j][2]; - MatrixAbsorb(a, b, c, dettri, matrix2); - - for (ii = 0; ii < 3; ii++) - { - for (jj = 0; jj < 3; jj++) - { left[ii, jj] = matrix1[ii, jj] + matrix2[ii, jj]; } - } - - for (ii = 0; ii < 3; ii++) - { - temp[ii] = 0; - for (k = 0; k < amesh.Ns; k++) - { temp[ii] += amesh.W[i][k][smesh.Region[smesh.T[j][ii]]] * flux[k][j][ii]; } - } - - source_corr = smesh.A[j] / 12; - SourceAssign(us[j], temp, right, RHS[i][j], dettri, source_corr); - - for (k = 0; k < 3; k++) - { - if (smesh.Bd2[bi][j][k] > 0) - { - for (ii = 0; ii < 2; ii++) - { - for (jj = 0; jj < 2; jj++) - { left[index[k, ii], index[k, jj]] += smesh.Bd2[bi][j][k] * bv[index[k, ii], index[k, jj]]; } - } - } - else if (smesh.Bd2[bi][j][k] < 0) - { - if (smesh.So2[j][k] > -1)// upwind flux from the internal reflection and boundary source - { - edge = smesh.So2[j][k]; - - temp[0] = 0; - temp[1] = 0; - for (m = 0; m < 2; m++) - { - for (ii = 0; ii < 2; ii++) - { - temp[ii] += flux[bb.Ri[edge][i][m]][j][smesh.E2[edge][ii]] * bb.Ri2[edge][i][m]; - } - for (ii = 0; ii < 2; ii++) - { - temp[ii] += q[bb.Si[edge][i][m]][edge][ii] * bb.Si2[edge][i][m]; - } - } - - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += temp[jj] * bv[smesh.E2[edge][ii], smesh.E2[edge][jj]]; - } - right[smesh.E2[edge][ii]] += -smesh.Bd2[bi][j][k] * tempd; - } - } - else // upwind flux from the adjacent triangle - { - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += flux[i][smesh.Bd[bi][j][3 * k]][smesh.Bd[bi][j][3 * k + jj + 1]] * bv[index[k, ii], index[k, jj]]; - } - right[index[k, ii]] += -smesh.Bd2[bi][j][k] * tempd; - } - } - } - } - // compute "right-left*flux" and then find the corresponding nodal values of "res". - tempd = 0; - for (ii = 0; ii < 3; ii++) - { - temp[ii] = (right[ii] - left[ii, 0] * flux[i][j][0] - left[ii, 1] * flux[i][j][1] - left[ii, 2] * flux[i][j][2]) / source_corr; - tempd += temp[ii]; - } - tempd = tempd / 4; - for (ii = 0; ii < 3; ii++) - { res[i][j][ii] = temp[ii] - tempd; } - } - } - } - else - { - for (i = 0; i < amesh.Ns; i++) - { - bi = i * aMeshLevel; - for (j = 0; j < smesh.Nt; j++) - { - dettri = 2 * smesh.A[j]; - cosi = amesh.Ang[i][0]; sini = amesh.Ang[i][1]; - a = cosi * (smesh.P[smesh.T[j][2]][1] - smesh.P[smesh.T[j][0]][1]) + sini * (smesh.P[smesh.T[j][0]][0] - smesh.P[smesh.T[j][2]][0]); - b = cosi * (smesh.P[smesh.T[j][0]][1] - smesh.P[smesh.T[j][1]][1]) + sini * (smesh.P[smesh.T[j][1]][0] - smesh.P[smesh.T[j][0]][0]); - MatrixConvec(a, b, matrix1); - - a = ua[j][0] + us[j][0]; - b = ua[j][1] + us[j][1]; - c = ua[j][2] + us[j][2]; - MatrixAbsorb(a, b, c, dettri, matrix2); - - for (ii = 0; ii < 3; ii++) - { - for (jj = 0; jj < 3; jj++) - { left[ii, jj] = matrix1[ii, jj] + matrix2[ii, jj]; } - } - - for (ii = 0; ii < 3; ii++) - { - temp[ii] = 0; - for (k = 0; k < amesh.Ns; k++) - { temp[ii] += amesh.W[i][k][smesh.Region[smesh.T[j][ii]]] * flux[k][j][ii]; } - } - - source_corr = smesh.A[j] / 12; - SourceAssign(us[j], temp, right, RHS[i][j], dettri, source_corr); - - - for (k = 0; k < 3; k++) - { - if (smesh.Bd2[bi][j][k] > 0) - { - for (ii = 0; ii < 2; ii++) - { - for (jj = 0; jj < 2; jj++) - { left[index[k, ii], index[k, jj]] += smesh.Bd2[bi][j][k] * bv[index[k, ii], index[k, jj]]; } - } - } - else if (smesh.Bd2[bi][j][k] < 0) - { - if (smesh.So2[j][k] > -1) - // "tri" is at the domain boundary: upwind flux is from boundary source only. - { - edge = smesh.So2[j][k]; - - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += q[i][edge][jj] * bv[smesh.E2[edge][ii], smesh.E2[edge][jj]]; - } - right[smesh.E2[edge][ii]] += -smesh.Bd2[bi][j][k] * tempd; - } - } - else - { - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += flux[i][smesh.Bd[bi][j][3 * k]][smesh.Bd[bi][j][3 * k + jj + 1]] * bv[index[k, ii], index[k, jj]]; - } - right[index[k, ii]] += -smesh.Bd2[bi][j][k] * tempd; - } - } - } - } - // compute "right-left*flux" and then find the corresponding nodal values of "res". - tempd = 0; - for (ii = 0; ii < 3; ii++) - { - temp[ii] = (right[ii] - left[ii, 0] * flux[i][j][0] - left[ii, 1] * flux[i][j][1] - left[ii, 2] * flux[i][j][2]) / source_corr; - tempd += temp[ii]; - } - tempd = tempd / 4; - for (ii = 0; ii < 3; ii++) - { res[i][j][ii] = temp[ii] - tempd; } - } - } - } - } - - - public void EdgeTri(int nt, double[] theta, double[][] p, int[][] p2, int[][] t, int[][] bd, double[][] bd2, int[][] so2) - // Purpose: this function is to compute edge integrals for each triangle. - // - // Let "n" be the normal of some edge of the triangle and "s" be the angular direction with - // the assumption for conforming mesh that one edge is shared by at most two triangles. - // - // Case 1 (s dot n>0): the edge flux is outgoing from the triangle. - // Case 2 (s dot n<0): - // Case 2.1 (internal triangle): the edge flux is incoming from the adjacent triangle. - // Case 2.2 (boundary triangle): the edge flux is incoming from the boundary source. - // - // For each angular direction "s" on each mesh "sMeshLevel", the edge integrals can be assembled from the following: - // - // bd[nt][9]: in Case 2.1, it saves the adjacent triangle of the current triangle, the local order of shared nodes - // in the adjacent triangle in the order of edges: index "0" to "2" for edge "1" in the current triangle, - // index "3" to "5" for edge "2" and index "6" to "8" for edge "3". - // Example: assume if the upwind flux for edge "2" of triangle "50" comes from triangle "10" with - // the local order "1" and "3", then bd[50][3]=10, bd[50][4]=1 and bd[50][5]=3. - // in Case 2.2, the corresponding boundary source is found through the function "boundary". - // - // bd2[nt][3]: the values of "s dot n * L (length of the edge)" by the order of edges of the current triangle. - // - // Note: we define "23" as the 1st edge, "31" as the 2nd and "12" as the 3rd. - { - int i, j; - double a = theta[0], b = theta[1]; // a=cos(theta), b=sin(theta) - double x1 = 0, y1 = 0, x2 = 0, y2 = 0, x3 = 0, y3 = 0, dx, dy, convConvConvTol = 1e-6; - for (i = 0; i < nt; i++) - { - x1 = p[t[i][0]][0]; y1 = p[t[i][0]][1]; - x2 = p[t[i][1]][0]; y2 = p[t[i][1]][1]; - x3 = p[t[i][2]][0]; y3 = p[t[i][2]][1]; - - // 1st edge "23" - dx = y3 - y2; dy = x2 - x3; - if (dx * (x1 - x3) + dy * (y1 - y3) > 0) // to make sure the normal of edge "23" is outgoing w.r.t. the triangle. - { dx = -dx; dy = -dy; } - - bd2[i][0] = a * dx + b * dy;// s dot n * L - if (Math.Abs(bd2[i][0]) < convConvConvTol) - { bd2[i][0] = 0; } - if (bd2[i][0] < 0 && so2[i][0] == -1) - // "bd2[i][n]<0 or s dot n <0" means that this edge has upwind flux from the other adjacent triangle or the boundary source. - // "so2[i][0]==-1" indicates that this triangle is not at the boundary. - { - bd[i][0] = FindTri(i, p2[t[i][1]], p2[t[i][2]]);// to find the adjacent triangle - if (bd[i][0] > -1) - { - for (j = 0; j < 3; j++) - { - if (t[bd[i][0]][j] == t[i][1])// to find the local order in the adjacent triangle of node "2" of the current triangle - { bd[i][1] = j; goto stop3; } - } - stop3: ; - for (j = 0; j < 3; j++) - { - if (t[bd[i][0]][j] == t[i][2])// to find the local order in the adjacent triangle of node "3" of the current triangle - { bd[i][2] = j; goto stop4; } - } - stop4: ; - } - } - - // 2nd edge "31" - dx = y1 - y3; dy = x3 - x1; - if (dx * (x2 - x3) + dy * (y2 - y3) > 0) - { dx = -dx; dy = -dy; } - - bd2[i][1] = a * dx + b * dy; - if (Math.Abs(bd2[i][1]) < convConvConvTol) - { bd2[i][1] = 0; } - if (bd2[i][1] < 0 && so2[i][1] == -1) - { - bd[i][3] = FindTri(i, p2[t[i][0]], p2[t[i][2]]); - if (bd[i][3] > -1) - { - for (j = 0; j < 3; j++) - { - if (t[bd[i][3]][j] == t[i][2]) - { bd[i][4] = j; goto stop5; } - } - stop5: ; - for (j = 0; j < 3; j++) - { - if (t[bd[i][3]][j] == t[i][0]) - { bd[i][5] = j; goto stop6; } - } - stop6: ; - } - } - - // 3rd edge "12" - dx = y2 - y1; dy = x1 - x2; - if (dx * (x3 - x1) + dy * (y3 - y1) > 0) - { dx = -dx; dy = -dy; } - - bd2[i][2] = a * dx + b * dy; - if (Math.Abs(bd2[i][2]) < convConvConvTol) - { bd2[i][2] = 0; } - if (bd2[i][2] < 0 && so2[i][2] == -1) - { - bd[i][6] = FindTri(i, p2[t[i][0]], p2[t[i][1]]); - if (bd[i][6] > -1) - { - for (j = 0; j < 3; j++) - { - if (t[bd[i][6]][j] == t[i][0]) - { bd[i][7] = j; goto stop1; } - } - stop1: ; - for (j = 0; j < 3; j++) - { - if (t[bd[i][6]][j] == t[i][1]) - { bd[i][8] = j; goto stop2; } - } - stop2: ; - } - } - } - } - - - - - - private int FindTri(int tri, int[] pt1, int[] pt2) - // Purpose: this function is find the adjacent triangle of the current triangle "tri" - // Note: the edge can not be at the boundary of the domain. - { - int i, j, tri2 = -1, flag = 0; - for (i = 1; i <= pt1[0]; i++) - { - for (j = 1; j <= pt2[0]; j++) - { - if (pt1[i] == pt2[j]) - { - if (pt1[i] != tri) - { - flag = 1; - tri2 = pt1[i]; - goto stop; - } - } - } - } - if (flag == 0) - { Console.WriteLine("nonconforming mesh!!!\n"); } - stop: return tri2; - } - - - private int FindTri2(int[] pt1, int[] pt2) - // Purpose: this function is to find the triangle containing the boundary. - // Note: "pt1" and "pt2" have to be boundary nodes of the same boundary. - { - int i, j, tri = -1; - for (i = 1; i <= pt1[0]; i++) - { - for (j = 1; j <= pt2[0]; j++) - { - if (pt1[i] == pt2[j]) - { - tri = pt1[i]; - goto stop; - } - } - } - stop: return tri; - } - - - public void FtoC(int nt_c, int ns_c, double[][][] flux, double[][][] cflux, int[][] smap, double[][][][] fc) - // Purpose: this function describes fine-to-coarse operator simultaneouly in angle and space, - // i.e., the combination of "ftoc_a" and "ftoc_s". - { - int i, j, k, m, n; - for (i = 0; i < ns_c; i++) - { - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - cflux[i][j][k] = 0; - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - cflux[i][j][k] += flux[2 * i][smap[j][m]][n] * fc[j][k][m - 1][n]; - } - } - } - } - } - } - - - public void FtoC_a(int nt, int ns_c, double[][][] flux, double[][][] cflux) - // Purpose: this function describes angular fine-to-coarse operator by simple restriction. - { - int i, j, k; - for (i = 0; i < nt; i++) - { - for (j = 0; j < ns_c; j++) - { - for (k = 0; k < 3; k++) - { - cflux[j][i][k] = flux[2 * j][i][k]; - - } - } - } - } - - - public void FtoC_s(int nt_c, int ns, double[][][] flux, double[][][] cflux, int[][] smap, double[][][][] fc) - // Purpose: this function describes spatial fine-to-coarse operator by L2 projection for piecewise linear DG. - { - int i, j, k, m, n; - for (i = 0; i < ns; i++) - { - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - cflux[i][j][k] = 0; - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - cflux[i][j][k] += flux[i][smap[j][m]][n] * fc[j][k][m - 1][n]; - - } - } - } - } - } - } - - - public void FtoC_s2(int nt_c, double[][] f, double[][] cf, int[][] smap, double[][][][] fc) - // Purpose: this function describes spatial fine-to-coarse operator by L2 projection for piecewise linear DG. - { - int j, k, m, n; - for (j = 0; j < nt_c; j++) - { - for (k = 0; k < 3; k++) - { - cf[j][k] = 0; - for (m = 1; m <= smap[j][0]; m++) - { - for (n = 0; n < 3; n++) - { - cf[j][k] += f[smap[j][m]][n] * fc[j][k][m - 1][n]; - } - } - } - } - } - - - private void MatrixAbsorb(double a, double b, double c, double dettri, double[,] matrix) - //Purpose of this matrix to assign absorption terms to the matrix - { - matrix[0, 0] = dettri / 60 * (3 * a + b + c); - matrix[0, 1] = dettri / 120 * (2 * a + 2 * b + c); - matrix[0, 2] = dettri / 120 * (2 * a + b + 2 * c); - matrix[1, 0] = matrix[0, 1]; - matrix[1, 1] = dettri / 60 * (a + 3 * b + c); - matrix[1, 2] = dettri / 120 * (a + 2 * b + 2 * c); - matrix[2, 0] = matrix[0, 2]; - matrix[2, 1] = matrix[1, 2]; - matrix[2, 2] = dettri / 60 * (a + b + 3 * c); - } - - - private void MatrixConvec(double a, double b, double[,] matrix) - //Purpose of this matrix to assign convection terms to the matrix - { - matrix[0, 0] = (a + b) / 6; - matrix[0, 1] = matrix[0, 0]; - matrix[0, 2] = matrix[0, 0]; - matrix[1, 0] = -a / 6; - matrix[1, 1] = matrix[1, 0]; - matrix[1, 2] = matrix[1, 0]; - matrix[2, 0] = -b / 6; - matrix[2, 1] = matrix[2, 0]; - matrix[2, 2] = matrix[2, 0]; - } - - - private void MatrixSolver(double[,] A, double[] B, double[] sol) - // Purpose: this function is to solve a 3 by 3 nonsingular linear system by Gaussian elimination with pivoting. - { - int[] ind = new int[3] { 0, 1, 2 }; - int i, j, temp2; - double temp; - - if (A[0, 0] == 0) - { - if (A[1, 0] == 0) - { - ind[0] = 2; ind[2] = 0; - } - else - { - ind[0] = 1; ind[1] = 0; - } - } - - temp = A[ind[0], 0]; - B[ind[0]] = B[ind[0]] / temp; - for (j = 0; j < 3; j++) - { A[ind[0], j] = A[ind[0], j] / temp; } - - for (i = 1; i < 3; i++) - { - if (A[ind[i], 0] == 0) - { } - else - { - temp = A[ind[i], 0]; - B[ind[i]] = B[ind[i]] / temp - B[ind[0]]; - for (j = 0; j < 3; j++) - { - A[ind[i], j] = A[ind[i], j] / temp - A[ind[0], j]; - } - } - } - - if (A[ind[1], 1] == 0.0) - { temp2 = ind[1]; ind[1] = ind[2]; ind[2] = temp2; } - temp = A[ind[1], 1]; - B[ind[1]] = B[ind[1]] / temp; - for (j = 1; j < 3; j++) - { A[ind[1], j] = A[ind[1], j] / temp; } - - if (A[ind[2], 1] == 0.0) - { } - else - { - temp = A[ind[2], 1]; - B[ind[2]] = B[ind[2]] / temp - B[ind[1]]; - for (j = 1; j < 3; j++) - { - A[ind[2], j] = A[ind[2], j] / temp - A[ind[1], j]; - } - } - - sol[2] = B[ind[2]] / A[ind[2], 2]; - sol[1] = (B[ind[1]] - sol[2] * A[ind[1], 2]) / A[ind[1], 1]; - sol[0] = (B[ind[0]] - sol[2] * A[ind[0], 2] - sol[1] * A[ind[0], 1]) / A[ind[0], 0]; - } - - - public double MgCycle(AngularMesh[] amesh, SpatialMesh[] smesh, BoundaryCoupling[] b, double[][][][] q, - double[][][][] RHS, double[][][] ua, double[][][] us, double[][][][] flux, double[][][][] d, - int n1, int n2, int aMeshLevel, int aMeshLevel0, int sMeshLevel, int sMeshLevel0, int NS, int vacuum, int mgMethod) - - // Purpose: this function contains the multigrid methods with V-cycle. - // AMG: angular multigrid method. - // SMG: spatial multigrid method. - // MG1: simultaneou angular and spatial multigrid method. - // MG2: angular multigrid first, then spatial multigrid. - // MG3: spatial multigrid first, then angular multigrid. - // MG4: alternative angular and spatial multigrid. - // MG4_a: angular multigrid first; - // MG4_s: spatial multigrid first. - { - int i, ns, nt, da, ds, level = -1; - double res = 1e10; - - nt = smesh[sMeshLevel].Nt; - ns = amesh[aMeshLevel].Ns; - ds = sMeshLevel - sMeshLevel0; - da = aMeshLevel - aMeshLevel0; - - switch (mgMethod) - { - case 1: //AMG: - level = da; - break; - case 2: //SMG: - level = ds; - break; - case 3: //MG1: - level = Math.Max(da, ds); - break; - case 4: //MG2: - level = ds + da; - break; - case 5: //MG3: - level = ds + da; - break; - case 6: //MG4_a: - level = ds + da; - break; - case 7: //MG4_s: - level = ds + da; - break; - } - - for (i = 0; i < n1; i++) - { - Relaxation(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], vacuum); - } - - switch (mgMethod) - { - case 1://AMG: - { - if (aMeshLevel == aMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - } - break; - case 2://SMG: - { - if (sMeshLevel == sMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - break; - case 3://MG1: - { - if (aMeshLevel == aMeshLevel0) - { - if (sMeshLevel == sMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - else - { - if (sMeshLevel == sMeshLevel0) - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC(smesh[sMeshLevel - 1].Nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF(smesh[sMeshLevel - 1].Nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - } - break; - case 4://MG2: - { - if (aMeshLevel == aMeshLevel0) - { - if (sMeshLevel == sMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - } - break; - case 5://MG3: - { - if (sMeshLevel == sMeshLevel0) - { - if (aMeshLevel == aMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - break; - case 6://MG4_a: - { - if (aMeshLevel == aMeshLevel0) - { - if (sMeshLevel == sMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - else - { - if (sMeshLevel == sMeshLevel0) - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - else - { - mgMethod = 7;//MG4_s - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - } - } - break; - case 7://MG4_s: - { - if (aMeshLevel == aMeshLevel0) - { - if (sMeshLevel == sMeshLevel0) - { } - else - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - else - { - if (sMeshLevel == sMeshLevel0) - { - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_a(nt, amesh[aMeshLevel - 1].Ns, d[level], RHS[level - 1]); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel - 1, aMeshLevel0, sMeshLevel, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_a(nt, amesh[aMeshLevel - 1].Ns, flux[level], flux[level - 1]); - } - else - { - mgMethod = 6; - Defect(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], d[level], vacuum); - res = Residual(nt, ns, d[level], smesh[sMeshLevel].A); - FtoC_s(smesh[sMeshLevel - 1].Nt, ns, d[level], RHS[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Fc); - MgCycle(amesh, smesh, b, q, RHS, ua, us, flux, d, n1, n2, aMeshLevel, aMeshLevel0, sMeshLevel - 1, sMeshLevel0, NS, vacuum, mgMethod); - CtoF_s(smesh[sMeshLevel - 1].Nt, ns, flux[level], flux[level - 1], smesh[sMeshLevel].Smap, smesh[sMeshLevel].Cf); - } - } - } - break; - } - - for (i = 0; i < n2; i++) - { - Relaxation(amesh[aMeshLevel], smesh[sMeshLevel], NS, RHS[level], ua[sMeshLevel], us[sMeshLevel], flux[level], b[level], q[level], vacuum); - } - - return res; - } - - private double Mod2pi(double x) - // Purpose: this function is to transfer angle "x" into [0 2*Pi) - { - double y; - if (x < 2 * Pi) - { - if (x < 0) - { - y = x + 2 * Pi; - } - else - { - y = x; - } - } - else - { - y = x - 2 * Pi; - } - return y; - } - - private double Reflection(double theta_i, double ni, double no) - // Purpose: this function is to find the reflection energy ratio for the reflected angle "theta_i" by tracing-back computation. - { - double r, theta_t, temp1, temp2; - if (Math.Abs(theta_i) < 1e-6) - { - temp1 = (ni - no) / (ni + no); - r = temp1 * temp1; - } - else - { - temp1 = Math.Sin(theta_i) * ni / no; - if (temp1 < 1) - { - theta_t = Math.Asin(temp1); - temp1 = Math.Sin(theta_i - theta_t) / Math.Sin(theta_i + theta_t); - temp2 = Math.Tan(theta_i - theta_t) / Math.Tan(theta_i + theta_t); - r = 0.5 * (temp1 * temp1 + temp2 * temp2); - } - else - { - r = 1; - } - } - return r; - } - - - private void Refraction(double[] temp, double theta_r, double ni, double no) - // Purpose: this function is to find the refraction energy ratio for the refracted angle "theta_r" by tracing-back computation. - { - double r, theta_i, temp1, temp2; - if (Math.Abs(theta_r) < 1e-6) - { - temp1 = (ni - no) / (ni + no); - r = 1 - temp1 * temp1; - theta_i = 0; - } - else - { - temp1 = Math.Sin(theta_r) * ni / no; - if (temp1 < 1) - { - theta_i = Math.Asin(temp1); - temp1 = Math.Sin(theta_i - theta_r) / Math.Sin(theta_i + theta_r); - temp2 = Math.Tan(theta_i - theta_r) / Math.Tan(theta_i + theta_r); - r = 1 - 0.5 * (temp1 * temp1 + temp2 * temp2); - } - else - { - theta_i = Pi / 2; - r = 0; - } - } - temp[0] = r; temp[1] = theta_i; - } - - private void Relaxation(AngularMesh amesh, SpatialMesh smesh, int Ns, double[][][] RHS, double[][] ua, double[][] us, double[][][] flux, - BoundaryCoupling bb, double[][][] q, int vacuum) - - // Purpose: this function is improved source-iteration (ISI) with vacuum or reflection boundary condition. - { - int i, j, k, m, ii, jj, tri, bi, edge = -1; - double[,] left = new double[3, 3]; - double[] right = new double[3]; - double[] temp = new double[3]; - double tempd; - double dettri, cosi, sini, a, b, c, sourceCorr; - double[,] bv = new double[3, 3] { { 1.0 / 3, 1.0 / 6, 1.0 / 6 }, { 1.0 / 6, 1.0 / 3, 1.0 / 6 }, { 1.0 / 6, 1.0 / 6, 1.0 / 3 } }; - double[,] matrix1 = new double[3, 3]; - double[,] matrix2 = new double[3, 3]; - - int[,] index = new int[3, 2]; - - int ns = amesh.Ns; - int aMeshLevel = Ns / ns; - int nt = smesh.Nt; - - index[0, 0] = 1; index[0, 1] = 2; - index[1, 0] = 2; index[1, 1] = 0; - index[2, 0] = 0; index[2, 1] = 1; - - - if (vacuum == 0)// reflection boundary condition - { - for (i = 0; i < ns; i++) - { - bi = i * aMeshLevel; - // "bi" is the angular index of the coarse angle on the fine angular mesh - // since sweep ordering is saved on the finest angular mesh for each spatial mesh for simplicity. - for (j = 0; j < nt; j++) - { - tri = smesh.So[bi][j];// the current sweep triangle by sweep ordering - dettri = 2 * smesh.A[tri]; - cosi = amesh.Ang[i][0]; sini = amesh.Ang[i][1]; - - // matrix1: convection term - a = cosi * (smesh.P[smesh.T[tri][2]][1] - smesh.P[smesh.T[tri][0]][1]) + sini * (smesh.P[smesh.T[tri][0]][0] - smesh.P[smesh.T[tri][2]][0]); - b = cosi * (smesh.P[smesh.T[tri][0]][1] - smesh.P[smesh.T[tri][1]][1]) + sini * (smesh.P[smesh.T[tri][1]][0] - smesh.P[smesh.T[tri][0]][0]); - MatrixConvec(a, b, matrix1); - - // matrix2: absorption term - a = ua[tri][0] + (1 - amesh.W[i][i][smesh.Region[smesh.T[j][0]]]) * us[tri][0]; - b = ua[tri][1] + (1 - amesh.W[i][i][smesh.Region[smesh.T[j][1]]]) * us[tri][1]; - c = ua[tri][2] + (1 - amesh.W[i][i][smesh.Region[smesh.T[j][2]]]) * us[tri][2]; - MatrixAbsorb(a, b, c, dettri, matrix2); - - // left[][]: convection+absorption - for (ii = 0; ii < 3; ii++) - { - for (jj = 0; jj < 3; jj++) - { left[ii, jj] = matrix1[ii, jj] + matrix2[ii, jj]; } - } - - // temp: scattering contribution to the direction "i" from all directions except "i". - for (ii = 0; ii < 3; ii++) - { - temp[ii] = 0; - for (k = 0; k < ns; k++) - { - temp[ii] += amesh.W[i][k][smesh.Region[smesh.T[j][ii]]] * flux[k][tri][ii]; - } - temp[ii] += -amesh.W[i][i][smesh.Region[smesh.T[j][ii]]] * flux[i][tri][ii]; - } - - sourceCorr = smesh.A[tri] / 12; - SourceAssign(us[tri], temp, right, RHS[i][tri], dettri, sourceCorr); - - - // add edge contributions to left, or add upwind fluxes to right from boundary source or the adjacent triangle - for (k = 0; k < 3; k++) - { // boundary outgoing flux (s dot n>0): add boundary contributions to left - if (smesh.Bd2[bi][tri][k] > 0) - { - for (ii = 0; ii < 2; ii++) - { - for (jj = 0; jj < 2; jj++) - { left[index[k, ii], index[k, jj]] += smesh.Bd2[bi][tri][k] * bv[index[k, ii], index[k, jj]]; } - } - } - // boundary incoming flux (s dot n<0): add upwind fluxes to right - else if (smesh.Bd2[bi][tri][k] < 0) - { - if (smesh.So2[tri][k] > -1) - // "tri" is at the domain boundary: upwind flux is from internal reflection and boundary source. - { - edge = smesh.So2[tri][k]; - - for (ii = 0; ii < 2; ii++) - { - temp[ii] = 0; - for (m = 0; m < 2; m++) - { - //temp[ii]+=flux[ri[edge][i][m]][tri][index[k][ii]]*ri2[edge][i][m]; - temp[ii] += flux[bb.Ri[edge][i][m]][tri][smesh.E2[edge][ii]] * bb.Ri2[edge][i][m]; - temp[ii] += q[bb.Si[edge][i][m]][edge][ii] * bb.Si2[edge][i][m]; - } - } - // note: we distinguish boundary source from internal source for reflection boundary condition - // due to refraction index mismatch at the boundary. - - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += temp[jj] * bv[smesh.E2[edge][ii], smesh.E2[edge][jj]]; - } - right[smesh.E2[edge][ii]] += -smesh.Bd2[bi][tri][k] * tempd; - } - } - else - // "tri" is not at the domain boundary: upwind flux is from the adjacent triangle. - { - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += flux[i][smesh.Bd[bi][tri][3 * k]][smesh.Bd[bi][tri][3 * k + jj + 1]] * bv[index[k, ii], index[k, jj]]; - } - right[index[k, ii]] += -smesh.Bd2[bi][tri][k] * tempd; - } - } - } - } - MatrixSolver(left, right, flux[i][tri]);// update the nodal values at "tri" - } - } - } - else // vacuum boundary condition - { - for (i = 0; i < ns; i++) - { - bi = i * aMeshLevel; - for (j = 0; j < nt; j++) - { - tri = smesh.So[bi][j]; - dettri = 2 * smesh.A[tri]; - cosi = amesh.Ang[i][0]; sini = amesh.Ang[i][1]; - a = cosi * (smesh.P[smesh.T[tri][2]][1] - smesh.P[smesh.T[tri][0]][1]) + sini * (smesh.P[smesh.T[tri][0]][0] - smesh.P[smesh.T[tri][2]][0]); - b = cosi * (smesh.P[smesh.T[tri][0]][1] - smesh.P[smesh.T[tri][1]][1]) + sini * (smesh.P[smesh.T[tri][1]][0] - smesh.P[smesh.T[tri][0]][0]); - MatrixConvec(a, b, matrix1); - - a = ua[tri][0] + (1 - amesh.W[i][i][smesh.Region[smesh.T[j][0]]]) * us[tri][0]; - b = ua[tri][1] + (1 - amesh.W[i][i][smesh.Region[smesh.T[j][1]]]) * us[tri][1]; - c = ua[tri][2] + (1 - amesh.W[i][i][smesh.Region[smesh.T[j][2]]]) * us[tri][2]; - MatrixAbsorb(a, b, c, dettri, matrix2); - - for (ii = 0; ii < 3; ii++) - { - for (jj = 0; jj < 3; jj++) - { left[ii, jj] = matrix1[ii, jj] + matrix2[ii, jj]; } - } - - for (ii = 0; ii < 3; ii++) - { - temp[ii] = 0; - for (k = 0; k < ns; k++) - { - temp[ii] += amesh.W[i][k][smesh.Region[smesh.T[j][ii]]] * flux[k][tri][ii]; - } - temp[ii] += -amesh.W[i][i][smesh.Region[smesh.T[j][ii]]] * flux[i][tri][ii]; - } - - sourceCorr = smesh.A[tri] / 12; - SourceAssign(us[tri], temp, right, RHS[i][tri], dettri, sourceCorr); - - - for (k = 0; k < 3; k++) - { - if (smesh.Bd2[bi][tri][k] > 0) - { - for (ii = 0; ii < 2; ii++) - { - for (jj = 0; jj < 2; jj++) - { left[index[k, ii], index[k, jj]] += smesh.Bd2[bi][tri][k] * bv[index[k, ii], index[k, jj]]; } - } - } - else if (smesh.Bd2[bi][tri][k] < 0) - { - if (smesh.So2[tri][k] > -1) - // "tri" is at the domain boundary: upwind flux is from boundary source only. - { - edge = smesh.So2[tri][k]; - - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += q[i][edge][jj] * bv[smesh.E2[edge][ii], smesh.E2[edge][jj]]; - } - right[smesh.E2[edge][ii]] += -smesh.Bd2[bi][tri][k] * tempd; - } - } - else - { - for (ii = 0; ii < 2; ii++) - { - tempd = 0; - for (jj = 0; jj < 2; jj++) - { - tempd += flux[i][smesh.Bd[bi][tri][3 * k]][smesh.Bd[bi][tri][3 * k + jj + 1]] * bv[index[k, ii], index[k, jj]]; - } - right[index[k, ii]] += -smesh.Bd2[bi][tri][k] * tempd; - } - } - } - } - MatrixSolver(left, right, flux[i][tri]); - } - } - } - } - - - public double Residual(int nt, int ns, double[][][] d, double[] a) - // Purpose: this function is to compute the residual. - { - double res = 0, temp; - int i, j, k; - - for (i = 0; i < nt; i++) - { - temp = 0; - for (j = 0; j < ns; j++) - { - for (k = 0; k < 3; k++) - { temp += Math.Abs(d[j][i][k]); } - } - res += temp * a[i]; - } - return res * 2 * Pi / ns / 3; - } - - - private void SourceAssign(double[] us, double[] temp, double[] right, double[] RHS, double dettri, double source_corr) - // Note: point (delta) source needs correction. - // right[][]: scattering+light source - { - double a, b, c, d, e, f; - - a = us[0]; - b = us[1]; - c = us[2]; - d = temp[0]; - e = temp[1]; - f = temp[2]; - - right[0] = dettri / 120 * (2 * a * (3 * d + e + f) + b * (2 * d + 2 * e + f) + c * (2 * d + e + 2 * f)) - + (2 * RHS[0] + RHS[1] + RHS[2]) * source_corr; - right[1] = dettri / 120 * (a * (2 * d + 2 * e + f) + 2 * b * (d + 3 * e + f) + c * (d + 2 * e + 2 * f)) - + (RHS[0] + 2 * RHS[1] + RHS[2]) * source_corr; - right[2] = dettri / 120 * (a * (2 * d + e + 2 * f) + b * (d + 2 * e + 2 * f) + 2 * c * (d + e + 3 * f)) - + (RHS[0] + RHS[1] + 2 * RHS[2]) * source_corr; - } - - - public void SpatialMapping(SpatialMesh cmesh, SpatialMesh fmesh, int[][] smap) - // Purpose: this function is to compute "cf" and "fc". - // - // cf[2][n][2]: coarse-to-fine spatial mapping of nodal values computed by linear interpolation - // fc[2][n][2]: fine-to-coarse spatial mapping of nodal values computed by L2 projection - // - // Note: 1D spatial mapping is straightforward due to the simple geometry, - // in particular, it is the same with respect to different intervals. - { - double[][] p; - double[][] c_f; - double[][] c_c; - double[] a; - double[] distance; - double area1, area2, area3, x, y, x1, y1, x2, y2, x3, y3; - int[][] t; - int nt_f, nt_c, i, j, tri, flag; - - p = cmesh.P; - t = cmesh.T; - c_c = cmesh.C; - c_f = fmesh.C; - a = cmesh.A; - nt_f = fmesh.Nt; - nt_c = cmesh.Nt; - - distance = new double[nt_c]; - - for (i = 0; i < nt_c; i++) - { smap[i][0] = 0; } - - for (i = 0; i < nt_f; i++) - { - flag = 0; - x = c_f[i][0]; y = c_f[i][1]; - for (j = 0; j < nt_c; j++) - { distance[j] = (x - c_c[j][0]) * (x - c_c[j][0]) + (y - c_c[j][1]) * (y - c_c[j][1]); } - for (j = 0; j < nt_c; j++) - { - tri = MathFunctions.FindMin(nt_c, distance);// find the triangle with minimal distance between (x,y) and centers of coarse triangles. - distance[tri] = 1e10; - x1 = p[t[tri][0]][0]; y1 = p[t[tri][0]][1]; - x2 = p[t[tri][1]][0]; y2 = p[t[tri][1]][1]; - x3 = p[t[tri][2]][0]; y3 = p[t[tri][2]][1]; - - area1 = MathFunctions.Area(x, y, x2, y2, x3, y3); - area2 = MathFunctions.Area(x1, y1, x, y, x3, y3); - area3 = MathFunctions.Area(x1, y1, x2, y2, x, y); - if (Math.Abs(area1 + area2 + area3 - a[tri]) / a[tri] < 1e-2)// If true, then the fine triangle "i" is in the coarse triangle "j". - { - flag = 1; - smap[tri][0] += 1; - smap[tri][smap[tri][0]] = i; - goto stop; - } - } - if (flag == 0) - { Console.WriteLine("spatial mapping has a problem\n"); } - stop: ; - } - } - - public void SpatialMapping2(SpatialMesh cmesh, SpatialMesh fmesh, int[][] smap, double[][][][] cf, double[][][][] fc, int tempsize) - - // Purpose: this function is to compute "cf" and "fc". - // - // cf[nt_c][3][smap[nt_c][0]][3]: coarse-to-fine spatial mapping of nodal values computed by linear interpolation - // fc[nt_c][3][smap[nt_c][0]][3]: fine-to-coarse spatial mapping of nodal values computed by L2 projection - // - // the 1st indexs the coarse triangles; - // the 2nd indexes 3 coarse nodes in the coarse triangle; - // the 3rd indexes the fine triangles contained in the coarse triangle; - // the 4th indexes 3 fine nodes in the fine triangle. - { - int[][] ind_p; - int np, nt_c, flag, i, j, k, m; - int[][] temp; - int[][] tf; - int[][] t; - double[][] temp2; - double tempd, x1, x2, x3, y1, y2, y3, x, y; - double[][] pf; - double[][] p; - double[] xf = new double[3]; - double[] yf = new double[3]; - double[] a; - double[] af; - - tf = fmesh.T; pf = fmesh.P; af = fmesh.A; - nt_c = cmesh.Nt; p = cmesh.P; t = cmesh.T; a = cmesh.A; - - temp = new int[3 * tempsize][]; - for (i = 0; i < 3 * tempsize; i++) - { - temp[i] = new int[2]; ; - for (j = 0; j < 2; j++) - { temp[i][j] = -1; } - } - - temp2 = new double[3 * tempsize][]; - for (i = 0; i < 3 * tempsize; i++) - { - temp2[i] = new double[3]; - for (j = 0; j < 3; j++) - { temp2[i][j] = 0; } - } - - ind_p = new int[tempsize][]; - for (j = 0; j < tempsize; j++) - { - ind_p[j] = new int[3]; - for (k = 0; k < 3; k++) - { ind_p[j][k] = -1; } - } - - for (i = 0; i < nt_c; i++) - { - // compute "ind_p" and "temp" - // temp[][]: find discrete (without repeat) nodes of the fine triangles contained in the coarse triangle - // and save it sequentially, i.e., order it in sequence as 1,2,3... and so on. - // the 1st indexes fine nodes; - // for the 2nd index, 1st entry is the global nodal index on fine mesh, 2nd entry is the repeatence of this node. - - // ind_p[][3]: the local nodal correspondence in fine triangles contained in the coarse triangle - // Example: if the 2nd node of the 3rd fine triangle is the 4th node as in "temp", - // then ind_p[3][2]=4. - - temp[0][0] = tf[smap[i][1]][0]; - temp[0][1] = 1; - np = 1; ind_p[0][0] = 0; - j = 0; - { - for (k = 1; k < 3; k++) - { - flag = 0;// flag indicates whether it is a new node that has not been found yet. - for (m = 0; m < np; m++) - { - if (tf[smap[i][j + 1]][k] == temp[m][0])// existing node - { - ind_p[j][k] = m; - temp[m][1] += 1; - flag = 1; - goto stop; - } - } - stop: ; - if (flag == 0)// new node - { - ind_p[j][k] = np; - temp[np][0] = tf[smap[i][j + 1]][k]; - temp[np][1] = 1; - np += 1; - } - } - } - - for (j = 1; j < smap[i][0]; j++) - { - for (k = 0; k < 3; k++) - { - flag = 0; - for (m = 0; m < np; m++) - { - if (tf[smap[i][j + 1]][k] == temp[m][0]) - { - ind_p[j][k] = m; - temp[m][1] += 1; - flag = 1; - goto stop2; - } - } - stop2: ; - if (flag == 0) - { - ind_p[j][k] = np; - temp[np][0] = tf[smap[i][j + 1]][k]; - temp[np][1] = 1; - np += 1; - } - } - } - - // compute temp2 - // temp2[][]: the nodal weights of non-repeat fine nodes in the coarse triangle. - x1 = p[t[i][0]][0]; y1 = p[t[i][0]][1]; - x2 = p[t[i][1]][0]; y2 = p[t[i][1]][1]; - x3 = p[t[i][2]][0]; y3 = p[t[i][2]][1]; - tempd = MathFunctions.Area(x1, y1, x2, y2, x3, y3); - for (j = 0; j < np; j++) - { - x = pf[temp[j][0]][0]; y = pf[temp[j][0]][1]; - temp2[j][1] = MathFunctions.Area(x1, y1, x, y, x3, y3) / tempd; - temp2[j][2] = MathFunctions.Area(x1, y1, x2, y2, x, y) / tempd; - temp2[j][0] = 1 - temp2[j][1] - temp2[j][2]; - } - - // compute "cf" - for (j = 0; j < smap[i][0]; j++) - { - for (k = 0; k < 3; k++) - { - for (m = 0; m < 3; m++) - { - cf[i][m][j][k] = temp2[ind_p[j][k]][m]; - } - } - } - - // compute "fc" - for (j = 0; j < smap[i][0]; j++) - { - for (k = 0; k < 3; k++) - { - xf[k] = temp2[ind_p[j][k]][1]; - yf[k] = temp2[ind_p[j][k]][2]; - } - - fc[i][1][j][0] = 2 * xf[0] + xf[1] + xf[2] - 1; - fc[i][2][j][0] = 2 * yf[0] + yf[1] + yf[2] - 1; - fc[i][0][j][0] = 1 - fc[i][1][j][0] - fc[i][2][j][0]; - - fc[i][1][j][1] = xf[0] + 2 * xf[1] + xf[2] - 1; - fc[i][2][j][1] = yf[0] + 2 * yf[1] + yf[2] - 1; - fc[i][0][j][1] = 1 - fc[i][1][j][1] - fc[i][2][j][1]; - - fc[i][1][j][2] = xf[0] + xf[1] + 2 * xf[2] - 1; - fc[i][2][j][2] = yf[0] + yf[1] + 2 * yf[2] - 1; - fc[i][0][j][2] = 1 - fc[i][1][j][2] - fc[i][2][j][2]; - - tempd = af[smap[i][j + 1]] / a[i]; - for (k = 0; k < 3; k++) - { - for (m = 0; m < 3; m++) - { fc[i][m][j][k] *= tempd; } - } - } - } - } - - } -} \ No newline at end of file diff --git a/src/Vts.FemModeling/MGRTE/2D/OutputCalculation.cs b/src/Vts.FemModeling/MGRTE/2D/OutputCalculation.cs deleted file mode 100644 index f025ae2fc..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/OutputCalculation.cs +++ /dev/null @@ -1,85 +0,0 @@ -using System; -using Vts.FemModeling.MGRTE._2D.DataStructures; - -namespace Vts.FemModeling.MGRTE._2D -{ - public class OutputCalculation - { - public double Pi = GlobalConstants.Pi; - - public Measurement RteOutput(double[][][] flux, double[][][] q, AngularMesh amesh, SpatialMesh smesh, BoundaryCoupling b, int vacuum) - - { - int i, j, k; - int nxy; - int[][] t; - int nt = smesh.Nt, ns = amesh.Ns, np = smesh.Np; - double dtheta = 2 * Pi / ns; - t = smesh.T; - - Measurement Det = new Measurement(); - - - nxy = (int)Math.Ceiling(Math.Sqrt(nt / 2.0)) + 1; - - Det.fluence = new double[np]; - - Det.radiance = new double[np][]; - for (i = 0; i < np; i++) - Det.radiance[i] = new double[ns]; - - Det.uxy = new double[nxy][]; - for (i = 0; i < nxy; i++) - Det.uxy[i] = new double[nxy]; - - Det.xloc = new double[nxy]; - Det.zloc = new double[nxy]; - Det.dx = new double[nxy]; - Det.dz = new double[nxy]; - Det.inten = new double[nxy*nxy]; - - - // compute radiance at each node - for (i = 0; i < nt; i++) - { - for (j = 0; j < ns; j++) - { - for (k = 0; k < 3; k++) - { - Det.radiance[t[i][k]][j] = flux[j][i][k]; - } - } - } - - // compute fluence at each node - for (i = 0; i < np; i++) - { - Det.fluence[i] = 0; - for (j = 0; j < ns; j++) - { - Det.fluence[i] += Det.radiance[i][j]; - } - Det.fluence[i] *= dtheta; - } - - MathFunctions.SquareTriMeshToGrid(ref smesh, ref Det.xloc, ref Det.zloc, ref Det.uxy, Det.fluence, nxy); - - for (i = 0; i < nxy; i++) - for (j = 0; j < nxy; j++) - { - Det.inten[i * nxy + j] = Det.uxy[i][j]; - } - - for (i = 0; i < nxy; i++) - { - Det.dx[i] = Det.xloc[1] - Det.xloc[0]; - Det.dz[i] = Det.zloc[1] - Det.zloc[0]; - } - - - - return Det; - } - } -} - diff --git a/src/Vts.FemModeling/MGRTE/2D/SolverMGRTE.cs b/src/Vts.FemModeling/MGRTE/2D/SolverMGRTE.cs deleted file mode 100644 index 88a717c71..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/SolverMGRTE.cs +++ /dev/null @@ -1,335 +0,0 @@ -using System; -using Vts.Common.Logging; -using Vts.FemModeling.MGRTE._2D.DataStructures; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Tissues; - -namespace Vts.FemModeling.MGRTE._2D -{ - /// - /// MG TRE Solver - /// - public class SolverMGRTE - { - /// - /// Execute MG RTE solver - /// - /// Simulation inputs - /// measurements - public static Measurement ExecuteMGRTE(SimulationInput input) - { - - int nMismatch; - int i, j, k, m, n; - int level; - double res = 0, res0 = 1, rho = 1.0; - int ds = input.MeshDataInput.SMeshLevel - input.SimulationOptionsInput.StartingSmeshLevel; - int da = input.MeshDataInput.AMeshLevel - input.SimulationOptionsInput.StartingAmeshLevel; - - - /* Read the initial time. */ - DateTime startTime1 = DateTime.Now; - ILogger logger = LoggerFactoryLocator.GetDefaultNLogFactory().Create(typeof (SolverMGRTE)); - - // step 1: compute "level" - // level: the indicator of mesh levels in multigrid - switch (input.SimulationOptionsInput.MethodMg) - { - case 1: - level = da; - break; - case 2: //SMG: - level = ds; - break; - case 3: //MG1: - level = Math.Max(da, ds); - break; - case 4: //MG2: - level = ds + da; - break; - case 5: //MG3: - level = ds + da; - break; - case 6: //MG4_a: - level = ds + da; - break; - case 7: //MG4_s: - level = ds + da; - break; - default: - level = -1; - break; - } - - - //Create Dynamic arrays based on above values - var amesh = new AngularMesh[input.MeshDataInput.AMeshLevel + 1]; - var smesh = new SpatialMesh[input.MeshDataInput.SMeshLevel + 1]; - var b = new BoundaryCoupling[level + 1]; - - var noflevel = new int[level + 1][]; - var ua = new double[input.MeshDataInput.SMeshLevel + 1][][]; - var us = new double[input.MeshDataInput.SMeshLevel + 1][][]; - var rhs = new double[level + 1][][][]; - var d = new double[level + 1][][][]; - var flux = new double[level + 1][][][]; - var q = new double[level + 1][][][]; - - var mgrid = new MultiGridCycle(); - var rteout = new OutputCalculation(); - - var tissueInput = (MultiEllipsoidTissueInput)input.TissueInput; - int incRegions = tissueInput.EllipsoidRegions.Length; - int tisRegions = tissueInput.LayerRegions.Length; - - double depth = 0.0; - for (i = 1; i < tissueInput.LayerRegions.Length - 1; i++) - depth += ((LayerTissueRegion)(tissueInput.LayerRegions[i])).ZRange.Stop - ((LayerTissueRegion)(tissueInput.LayerRegions[i])).ZRange.Start; - - input.MeshDataInput.SideLength = depth; - int totRegions = incRegions + tisRegions ; - - // MG-RTE does not converge when g = 1.0; - for (i = 0; i < totRegions; i++) - if (input.TissueInput.Regions[i].RegionOP.G >= 1.0) - input.TissueInput.Regions[i].RegionOP.G = 1.0 - 1e-5; - - // Check refractive index mismatch - if (Math.Abs(input.TissueInput.Regions[0].RegionOP.N - input.TissueInput.Regions[1].RegionOP.N) - / input.TissueInput.Regions[0].RegionOP.N < 0.01) // refraction index mismatch at the boundary - nMismatch = 1; - else - nMismatch = 0; - - - //Create spatial and angular mesh - MathFunctions.CreateSquareMesh(ref smesh, input.MeshDataInput.SMeshLevel, depth); - MathFunctions.AssignRegions(ref smesh, input.MeshDataInput.SMeshLevel, tissueInput); - MathFunctions.CreateAnglularMesh(ref amesh, input.MeshDataInput.AMeshLevel, tissueInput); - - - MathFunctions.SweepOrdering(ref smesh, amesh, input.MeshDataInput.SMeshLevel, input.MeshDataInput.AMeshLevel); - MathFunctions.SetMus(ref us, smesh, input); - MathFunctions.SetMua(ref ua, smesh, input); - - // load optical property, angular mesh, and spatial mesh files - Initialization.Initial( - ref amesh, ref smesh, ref flux, ref d, - ref rhs, ref q, ref noflevel, ref b, - level, input.SimulationOptionsInput.MethodMg,nMismatch,input.SimulationOptionsInput.NExternal, - input.SimulationOptionsInput.NExternal,input.MeshDataInput.AMeshLevel, input.SimulationOptionsInput.StartingAmeshLevel, - input.MeshDataInput.SMeshLevel, input.SimulationOptionsInput.StartingSmeshLevel, ua, us, mgrid); - - //Assign external source if available - if (input.ExtSourceInput != null) - { - IExtSource extsource = FemSourceFactory.GetExtSource(input.ExtSourceInput); - extsource.AssignMeshForExtSource(amesh, input.MeshDataInput.AMeshLevel, smesh, input.MeshDataInput.SMeshLevel, level, q); - } - - //Assign internal source if available - if (input.IntSourceInput != null) - { - //Assign an internal source - IIntSource intsource = FemSourceFactory.GetIntSource(input.IntSourceInput); - intsource.AssignMeshForIntSource(amesh, input.MeshDataInput.AMeshLevel, smesh, input.MeshDataInput.SMeshLevel, level, rhs); - } - - /* Read the end time. */ - DateTime stopTime1 = DateTime.Now; - /* Compute and print the duration of this first task. */ - TimeSpan duration1 = stopTime1 - startTime1; - - logger.Info(() => "Initlalization for RTE_2D takes " + duration1.TotalSeconds + " seconds\n"); - - //step 2: RTE solver - DateTime startTime2 = DateTime.Now; - - int ns = amesh[input.MeshDataInput.AMeshLevel].Ns; - - - if (input.SimulationOptionsInput.FullMg == 1) - { - int nt1, ns1; - int nt2 = smesh[noflevel[level][0]].Nt; - int ns2 = amesh[noflevel[level][1]].Ns; - - for (n = level - 1; n >= 0; n--) - { - nt1 = smesh[noflevel[n][0]].Nt; - ns1 = amesh[noflevel[n][1]].Ns; - if (nt1 == nt2) - { - mgrid.FtoC_a(nt1, ns1, rhs[n + 1], rhs[n]); - } - else - { - if (ns1 == ns2) - { - mgrid.FtoC_s(nt1, ns1, rhs[n + 1], rhs[n], smesh[noflevel[n][0] + 1].Smap, smesh[noflevel[n][0] + 1].Fc); - } - else - { - mgrid.FtoC(nt1, ns1, rhs[n + 1], rhs[n], smesh[noflevel[n][0] + 1].Smap, smesh[noflevel[n][0] + 1].Fc); - } - - } - nt2 = nt1; ns2 = ns1; - } - - nt1 = smesh[noflevel[0][0]].Nt; - ns1 = amesh[noflevel[0][1]].Ns; - - for (n = 0; n < level; n++) - { - if (input.SimulationOptionsInput.MethodMg == 6) - { - if (((level - n) % 2) == 0) - { - for (i = 0; i < input.SimulationOptionsInput.NCycle; i++) - { - res = mgrid.MgCycle(amesh, smesh, b, q, rhs, ua, us, flux, d, input.SimulationOptionsInput.NPreIterations, input.SimulationOptionsInput.NPostIterations, - noflevel[n][1], input.SimulationOptionsInput.StartingAmeshLevel, noflevel[n][0], input.SimulationOptionsInput.StartingSmeshLevel, ns, nMismatch, 6); - } - } - else - { - for (i = 0; i < input.SimulationOptionsInput.NCycle; i++) - { - mgrid.MgCycle(amesh, smesh, b, q, rhs, ua, us, flux, d, input.SimulationOptionsInput.NPreIterations, input.SimulationOptionsInput.NPostIterations, - noflevel[n][1], input.SimulationOptionsInput.StartingAmeshLevel, noflevel[n][0], input.SimulationOptionsInput.StartingSmeshLevel, ns, nMismatch, 7); - } - } - } - else - { - if (input.SimulationOptionsInput.MethodMg== 7) - { - if (((level - n) % 2) == 0) - { - for (i = 0; i < input.SimulationOptionsInput.NCycle; i++) - { - mgrid.MgCycle(amesh, smesh, b, q, rhs, ua, us, flux, d, input.SimulationOptionsInput.NPreIterations, input.SimulationOptionsInput.NPostIterations, - noflevel[n][1], input.SimulationOptionsInput.StartingAmeshLevel, noflevel[n][0], input.SimulationOptionsInput.StartingSmeshLevel, ns, nMismatch, 7); - } - } - else - { - for (i = 0; i < input.SimulationOptionsInput.NCycle; i++) - { - mgrid.MgCycle(amesh, smesh, b, q, rhs, ua, us, flux, d, input.SimulationOptionsInput.NPreIterations, input.SimulationOptionsInput.NPostIterations, - noflevel[n][1], input.SimulationOptionsInput.StartingAmeshLevel, noflevel[n][0], input.SimulationOptionsInput.StartingSmeshLevel, ns, nMismatch, 6); - } - } - } - else - { - for (i = 0; i < input.SimulationOptionsInput.NCycle; i++) - { - mgrid.MgCycle(amesh, smesh, b, q, rhs, ua, us, flux, d, input.SimulationOptionsInput.NPreIterations, input.SimulationOptionsInput.NPostIterations, - noflevel[n][1], input.SimulationOptionsInput.StartingAmeshLevel, noflevel[n][0], input.SimulationOptionsInput.StartingSmeshLevel, ns, nMismatch, - input.SimulationOptionsInput.MethodMg); - } - } - } - - nt2 = smesh[noflevel[n + 1][0]].Nt; - ns2 = amesh[noflevel[n + 1][1]].Ns; - if (nt1 == nt2) - { - mgrid.CtoF_a(nt1, ns1, flux[n + 1], flux[n]); - } - else - { - if (ns1 == ns2) - { - mgrid.CtoF_s(nt1, ns1, flux[n + 1], flux[n], smesh[noflevel[n][0] + 1].Smap, smesh[noflevel[n][0] + 1].Cf); - } - else - { - mgrid.CtoF(nt1, ns1, flux[n + 1], flux[n], smesh[noflevel[n][0] + 1].Smap, smesh[noflevel[n][0] + 1].Cf); - } - } - nt1 = nt2; ns1 = ns2; - for (m = 0; m <= n; m++) - { - for (i = 0; i < amesh[noflevel[m][1]].Ns; i++) - { - for (j = 0; j < smesh[noflevel[m][0]].Nt; j++) - { - for (k = 0; k < 3; k++) - { - flux[m][i][j][k] = 0; - } - } - } - } - } - } - - // 2.2. multigrid solver on the finest mesh - n = 0; - - while (n < input.SimulationOptionsInput.NIterations) - { - n++; - res = mgrid.MgCycle(amesh, smesh, b, q, rhs, ua, us, flux, d, input.SimulationOptionsInput.NPreIterations, input.SimulationOptionsInput.NPostIterations, input.MeshDataInput.AMeshLevel, - input.SimulationOptionsInput.StartingAmeshLevel, input.MeshDataInput.SMeshLevel, input.SimulationOptionsInput.StartingSmeshLevel, ns, nMismatch, input.SimulationOptionsInput.MethodMg); - for (m = 0; m < level; m++) - { - for (i = 0; i < amesh[noflevel[m][1]].Ns; i++) - { - for (j = 0; j < smesh[noflevel[m][0]].Nt; j++) - { - for (k = 0; k < 3; k++) - { - flux[m][i][j][k] = 0; - } - } - } - } - - - if (n > 1) - { - rho *= res / res0; - logger.Info(() => "Iteration: " + n + ", Current tolerance: " + res + "\n"); - - if (res < input.SimulationOptionsInput.ConvTolerance) - { - rho = Math.Pow(rho, 1.0 / (n - 1)); - n = input.SimulationOptionsInput.NIterations; - } - } - else - { - logger.Info(() => "Iteration: " + n + ", Current tolerance: " + res + "\n"); - res0 = res; - if (res < input.SimulationOptionsInput.ConvTolerance) - n = input.SimulationOptionsInput.NIterations; - } - - } - - // 2.3. compute the residual - //Mgrid.Defect(amesh[para.AMeshLevel], smesh[para.SMeshLevel], ns, RHS[level], ua[para.SMeshLevel], us[para.SMeshLevel], - // flux[level], b[level], q[level], d[level], vacuum); - //res = Mgrid.Residual(smesh[para.SMeshLevel].nt, amesh[para.AMeshLevel].ns, d[level], smesh[para.SMeshLevel].a); - - /* Read the start time. */ - DateTime stopTime2 = DateTime.Now; - TimeSpan duration2 = stopTime2 - startTime2; - TimeSpan duration3 = stopTime2 - startTime1; - - logger.Info(() => "Iteration time: " + duration2.TotalSeconds + "seconds\n"); - logger.Info(() => "Total time: " + duration3.TotalSeconds + "seconds, Final residual: " + res + "\n"); - - // step 3: postprocessing - // 3.1. output - Measurement measurement = rteout.RteOutput(flux[level], q[level], amesh[input.MeshDataInput.AMeshLevel], smesh[input.MeshDataInput.SMeshLevel], b[level], nMismatch); - - return measurement; - } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Sources/ExternalSources/ExtLineSource.cs b/src/Vts.FemModeling/MGRTE/2D/Sources/ExternalSources/ExtLineSource.cs deleted file mode 100644 index d0364671e..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Sources/ExternalSources/ExtLineSource.cs +++ /dev/null @@ -1,57 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// External Line source - /// - public class ExtLineSource : IExtSource - { - /// - /// General constructor for external line source - /// - /// - /// - /// - public ExtLineSource( - DoubleRange start, - DoubleRange end, - DoubleRange thetaRange) - { - Start = start; - End = end; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for external line source - /// - public ExtLineSource() - : this( - new DoubleRange(-0.25, 0), - new DoubleRange(0.25, 0), - new DoubleRange(-0.5 * Math.PI, 0.5 * Math.PI)) { } - - /// - /// Starting coordinates (x,z) of the line - /// - public DoubleRange Start { get; set; } - - /// - /// Ending coordinates (x,z) of the line - /// - public DoubleRange End { get; set; } - - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - - public void AssignMeshForExtSource(AngularMesh[] amesh, int ameshLevel, SpatialMesh[] smesh, int smeshLevel, int level, double[][][][] q) - { - } - } - -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Sources/ExternalSources/ExtPointSource.cs b/src/Vts.FemModeling/MGRTE/2D/Sources/ExternalSources/ExtPointSource.cs deleted file mode 100644 index 415b47968..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Sources/ExternalSources/ExtPointSource.cs +++ /dev/null @@ -1,106 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// External point source - /// - public class ExtPointSource : IExtSource - { - /// - /// General constructor for external point source - /// - /// - /// - public ExtPointSource( - DoubleRange launchPoint, - DoubleRange thetaRange) - { - LaunchPoint = launchPoint; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for external point source - /// - public ExtPointSource() - : this( - new DoubleRange(0, 0), - new DoubleRange(Math.PI, 2.0 * Math.PI)) { } - - /// - /// Launching coordinates (x,z) of the boundary point source - /// - public DoubleRange LaunchPoint { get; set; } - - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - /// - /// - /// - /// - /// - /// - /// - /// - /// - public void AssignMeshForExtSource(AngularMesh[] aMesh, int aMeshLevel, SpatialMesh[] sMesh, int sMeshLevel, int level, double[][][][] q) - { - double[] distance = new double[sMesh[sMeshLevel].Nt]; - int i, k; - double x,x1,x2,z,z1,z2; - int ns_start, ns_stop; - - x = LaunchPoint.Start; - z = LaunchPoint.Stop; - - ns_start = (int)(0.5 * ThetaRange.Start * aMesh[aMeshLevel].Ns/ Math.PI); - ns_stop = (int)(0.5 * ThetaRange.Stop * aMesh[aMeshLevel].Ns / Math.PI); - - for (i = 0; i < sMesh[sMeshLevel].Ne; i++) - { - x1 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].E[i][1]][0]; - x2 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].E[i][2]][0]; - z1 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].E[i][1]][1]; - z2 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].E[i][2]][1]; - - if ((z1 == z) && (z2 == z)) - { - if ((x1 - x >= 0) && (x - x2 >= 0)) - { - for (k = 0; k < aMesh[aMeshLevel].Ns; k++) - { - if ((k < ns_start) || (k > ns_stop)) - q[level][aMesh[aMeshLevel].Ns-1-k][i][0] = 0; - else - { - q[level][aMesh[aMeshLevel].Ns - 1 - k][i][0] = 0.5*Math.Abs((x2 - x)/(x1 - x2)); - q[level][aMesh[aMeshLevel].Ns - 1 - k][i][1] = 0.5*Math.Abs((x1 - x)/(x1 - x2)); - - } - } - } - if ((x2 - x >= 0) && (x - x1 >= 0)) - { - for (k = 0; k < aMesh[aMeshLevel].Ns; k++) - { - if ((k < ns_start) || (k > ns_stop)) - q[level][aMesh[aMeshLevel].Ns - 1 - k][i][0] = 0; - else - { - q[level][aMesh[aMeshLevel].Ns - 1 - k][i][0] = 0.5 * Math.Abs((x2 - x) / (x1 - x2)); - q[level][aMesh[aMeshLevel].Ns - 1 - k][i][1] = 0.5 * Math.Abs((x1 - x) / (x1 - x2)); - - } - } - } - } - } - } - } - -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DCircularSource.cs b/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DCircularSource.cs deleted file mode 100644 index 0cf21b57b..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DCircularSource.cs +++ /dev/null @@ -1,56 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Internal 2D Circular source - /// - public class Int2DCircularSource : IIntSource - { - /// - /// General constructor for 2D circular source - /// - /// - /// - /// - public Int2DCircularSource( - double radius, - DoubleRange center, - DoubleRange thetaRange) - { - Radius = radius; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D circular source - /// - public Int2DCircularSource() - : this( - 0.5, - new DoubleRange(0.5, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Radius of the circular source - /// - public double Radius { get; set; } - - /// - /// Center cooordinates (x,z) of the geometry - /// - public DoubleRange Center { get; set; } - - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - public void AssignMeshForIntSource(AngularMesh[] amesh, int ameshLevel, SpatialMesh[] smesh, int smeshLevel, int level, double[][][][] rhs) - { - } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DEllipticalSource.cs b/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DEllipticalSource.cs deleted file mode 100644 index 62ffc064f..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DEllipticalSource.cs +++ /dev/null @@ -1,71 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Internal 2D Elliptical source - /// - public class Int2DEllipticalSource : IIntSource - { - /// - /// General constructor for 2D elliptical source - /// - /// - /// - /// - /// - public Int2DEllipticalSource( - double aParameter, - double bParameter, - DoubleRange center, - DoubleRange thetaRange) - { - AParameter = aParameter; - BParameter = bParameter; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D elliptical source - /// - public Int2DEllipticalSource() - : this( - 0.5, - 0.25, - new DoubleRange(0.5, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// a Parameter of Elliptical source - /// - public double AParameter { get; set; } - /// - /// b Parameter of Elliptical source - /// - public double BParameter { get; set; } - /// - /// Center cooordinates (x,z) of the geometry - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - /// - /// - /// - /// - /// - /// - /// - /// - /// - public void AssignMeshForIntSource(AngularMesh[] amesh, int ameshLevel, SpatialMesh[] smesh, int smeshLevel, int level, double[][][][] rhs) - { - } - } - -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DPointSource.cs b/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DPointSource.cs deleted file mode 100644 index 84361c5a5..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DPointSource.cs +++ /dev/null @@ -1,121 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Internal 2D Point source - /// - public class Int2DPointSource : IIntSource - { - /// - /// General constructor for 2D point source - /// - /// - /// - public Int2DPointSource( - DoubleRange center, - DoubleRange thetaRange) - { - Center = center; - ThetaRange = thetaRange; - } - - /// - /// Default constructor for 2D point source - /// - public Int2DPointSource() - : this( - new DoubleRange(0, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Center cooordinates (x,z) of the geometry - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - /// - /// - /// - /// - /// - /// - /// - /// - /// - public void AssignMeshForIntSource(AngularMesh[] aMesh, int aMeshLevel, SpatialMesh[] sMesh, int sMeshLevel, int level, double[][][][] rhs) - { - double[] distance = new double[sMesh[sMeshLevel].Nt]; - int i,j,k; - double x,x1,x2,x3,z,z1,z2,z3; - int ns_start, ns_stop; - - - double[] area = new double[3]; - double area_total; - - x = Center.Start; - z = Center.Stop; - - ns_start = (int)(0.5*ThetaRange.Start * aMesh[aMeshLevel].Ns/ Math.PI); - ns_stop = (int)(0.5 * ThetaRange.Stop * aMesh[aMeshLevel].Ns / Math.PI); - - for (i = 0; i < sMesh[sMeshLevel].Nt; i++) - { - x1 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].T[i][0]][0]; - x2 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].T[i][1]][0]; - x3 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].T[i][2]][0]; - z1 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].T[i][0]][1]; - z2 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].T[i][1]][1]; - z3 = sMesh[sMeshLevel].P[sMesh[sMeshLevel].T[i][2]][1]; - - double detT = (z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3); - - double l1 = ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) / detT; - double l2 = ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) / detT; - double l3 = 1 - l1 - l2; - - - if (l1 >= 0.0) - { - if (l1 <= 1.0) - { - if (l2 >= 0.0) - { - if (l2 <= 1.0) - { - if (l3 >= 0.0) - { - if (l3 <= 1.0) - { - area[0] = MathFunctions.Area(x, z, x2, z2, x3, z3); - area[1] = MathFunctions.Area(x1, z1, x, z, x3, z3); - area[2] = MathFunctions.Area(x1, z1, x2, z2, x, z); - area_total = area[0] + area[1] + area[2]; - for (j = 0; j < 3; j++) - { - for (k = 0; k < aMesh[aMeshLevel].Ns; k++) - { - if ((k < ns_start) || (k > ns_stop)) - rhs[level][aMesh[aMeshLevel].Ns-1-k][i][j] = 0; - else - rhs[level][aMesh[aMeshLevel].Ns-1-k][i][j] = area[j] / area_total; - } - } - } - } - } - } - } - } - - - } - } - - } -} diff --git a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DRectangularSource.cs b/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DRectangularSource.cs deleted file mode 100644 index 1d417385e..000000000 --- a/src/Vts.FemModeling/MGRTE/2D/Sources/InternalSources/Int2DRectangularSource.cs +++ /dev/null @@ -1,63 +0,0 @@ -using System; -using Vts.Common; - -namespace Vts.FemModeling.MGRTE._2D.DataStructures -{ - /// - /// Internal 2D Rectangular source - /// - public class Int2DRectangularSource : IIntSource - { - /// - /// - /// - /// - /// - /// - /// - public Int2DRectangularSource( - double xLength, - double zHeight, - DoubleRange center, - DoubleRange thetaRange) - { - XLength = xLength; - ZHeight = zHeight; - Center = center; - ThetaRange = thetaRange; - } - - /// - /// - /// - public Int2DRectangularSource() - : this( - 0.5, - 0.5, - new DoubleRange(0.5, 0.5), - new DoubleRange(0, 2 * Math.PI)) { } - - /// - /// Length of the Rectangular source - /// - public double XLength { get; set; } - /// - /// Height of the Rectangular source - /// - public double ZHeight { get; set; } - /// - /// Center cooordinates (x,z) of the geometry - /// - public DoubleRange Center { get; set; } - /// - /// Theta angle range - /// - public DoubleRange ThetaRange { get; set; } - - - public void AssignMeshForIntSource(AngularMesh[] amesh, int ameshLevel, SpatialMesh[] smesh, int smeshLevel, int level, double[][][][] rhs) - { - } - - } -} diff --git a/src/Vts.FemModeling/Vts.FemModeling.csproj b/src/Vts.FemModeling/Vts.FemModeling.csproj deleted file mode 100644 index 2752294ec..000000000 --- a/src/Vts.FemModeling/Vts.FemModeling.csproj +++ /dev/null @@ -1,11 +0,0 @@ - - - - net6.0 - - - - - - - diff --git a/src/Vts.Mgrte.Console/NLog.config b/src/Vts.Mgrte.Console/NLog.config deleted file mode 100644 index a51a99bf9..000000000 --- a/src/Vts.Mgrte.Console/NLog.config +++ /dev/null @@ -1,14 +0,0 @@ - - - - - - - - - - - - - \ No newline at end of file diff --git a/src/Vts.Mgrte.Console/Program.cs b/src/Vts.Mgrte.Console/Program.cs deleted file mode 100644 index 6b8754903..000000000 --- a/src/Vts.Mgrte.Console/Program.cs +++ /dev/null @@ -1,16 +0,0 @@ -using Vts.FemModeling.MGRTE._2D; - -namespace Vts.Mgrte.Console -{ - class Program - { - static void Main(string[] args) - { - var input = new SimulationInput(); - - SolverMGRTE.ExecuteMGRTE(input); - - System.Console.ReadLine(); - } - } -} diff --git a/src/Vts.Mgrte.Console/Vts.Mgrte.Console.csproj b/src/Vts.Mgrte.Console/Vts.Mgrte.Console.csproj deleted file mode 100644 index cc9055d16..000000000 --- a/src/Vts.Mgrte.Console/Vts.Mgrte.Console.csproj +++ /dev/null @@ -1,20 +0,0 @@ - - - - Exe - net6.0 - mgrte - Vts.Mgrte.Console - - - - - - - - - Always - - - - diff --git a/src/Vts.sln b/src/Vts.sln index c1a3609eb..096439421 100644 --- a/src/Vts.sln +++ b/src/Vts.sln @@ -23,10 +23,6 @@ Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.MonteCarlo.PostProcesso EndProject Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Benchmark", "Vts.Benchmark\Vts.Benchmark.csproj", "{3B0676C4-5CF1-4A0B-97E4-0CDF2B035702}" EndProject -Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Mgrte.Console", "Vts.Mgrte.Console\Vts.Mgrte.Console.csproj", "{197229DB-8CEC-48DD-A953-FBC58D1A3EBC}" -EndProject -Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.FemModeling", "Vts.FemModeling\Vts.FemModeling.csproj", "{6BA37669-7BA8-4FF2-911A-0E70EECD450B}" -EndProject Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Scripting", "Vts.Scripting\Vts.Scripting.csproj", "{2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}" EndProject Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Scripting.Test", "Vts.Scripting.Test\Vts.Scripting.Test.csproj", "{529AE3EE-28AF-4C0F-8D07-DC89833128DA}" @@ -225,70 +221,6 @@ Global {3B0676C4-5CF1-4A0B-97E4-0CDF2B035702}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU {3B0676C4-5CF1-4A0B-97E4-0CDF2B035702}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU {3B0676C4-5CF1-4A0B-97E4-0CDF2B035702}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|Any CPU.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|Any CPU.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|Mixed Platforms.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|Mixed Platforms.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|Win32.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|Win32.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|x86.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Benchmark|x86.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|Any CPU.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|Any CPU.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|Win32.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|Win32.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|x86.ActiveCfg = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Debug|x86.Build.0 = Debug|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|Any CPU.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|Any CPU.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|Mixed Platforms.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|Win32.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|Win32.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|x86.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.Release|x86.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|Any CPU.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|Any CPU.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|Mixed Platforms.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|Mixed Platforms.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|Win32.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|Any CPU.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|Any CPU.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|Mixed Platforms.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|Mixed Platforms.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|Win32.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|Win32.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|x86.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Benchmark|x86.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|Any CPU.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|Any CPU.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|Win32.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|Win32.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|x86.ActiveCfg = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Debug|x86.Build.0 = Debug|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|Any CPU.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|Any CPU.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|Mixed Platforms.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|Win32.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|Win32.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|x86.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.Release|x86.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Any CPU.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Any CPU.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Mixed Platforms.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Mixed Platforms.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Win32.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU - {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Any CPU.ActiveCfg = Debug|Any CPU {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Any CPU.Build.0 = Debug|Any CPU {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Mixed Platforms.ActiveCfg = Debug|Any CPU @@ -365,13 +297,11 @@ Global {36DE9584-F033-4EC0-A18F-1AE60EE55754} = {920000AB-D955-4399-AD00-78523C740E2A} {8A1194A0-C773-483B-BD78-9D013105CEA1} = {920000AB-D955-4399-AD00-78523C740E2A} {3B0676C4-5CF1-4A0B-97E4-0CDF2B035702} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E} - {197229DB-8CEC-48DD-A953-FBC58D1A3EBC} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E} - {6BA37669-7BA8-4FF2-911A-0E70EECD450B} = {E87433B1-D691-4D36-BB6B-3B766192CC8E} {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E} {529AE3EE-28AF-4C0F-8D07-DC89833128DA} = {920000AB-D955-4399-AD00-78523C740E2A} EndGlobalSection GlobalSection(ExtensibilityGlobals) = postSolution - SolutionGuid = {BB1B236A-1BE9-476A-A4B3-8C0F51F1FDA7} EnterpriseLibraryConfigurationToolBinariesPath = packages\Unity.2.1.505.2\lib\NET35 + SolutionGuid = {BB1B236A-1BE9-476A-A4B3-8C0F51F1FDA7} EndGlobalSection EndGlobal diff --git a/src/packages/repositories.config b/src/packages/repositories.config deleted file mode 100644 index a4497b4a3..000000000 --- a/src/packages/repositories.config +++ /dev/null @@ -1,15 +0,0 @@ - - - - - - - - - - - - - - - \ No newline at end of file