Hydrodynamics Kernel (12Apr Version)
Jeffrey Oldham
oldham at codesourcery.com
Thu Apr 12 23:10:56 UTC 2001
Yesterday, I delivered a hydrodynamics kernel for testing the Field
class and to use for timings. Today, I replaced the two sequential
computations with parallel versions:
1) changing the initialization of the coordinates Fields
2) changing the computation of corner forces.
There are now no significant sequential computations in the code.
(There are a few; just count the number of semicolons. :)
Attached is
1) the main hydrodynamics kernel file (also included in the archive)
2) a uuencoded gzipped tar archive with all the source files
3) a patch with the differences between yesterday and today's source code
To compile the program,
a. uudecode foo.uu # produces hydrodynamics.tgz
b. tar xzvf hydrodynamics.tgz # extracts the five files
c. In Makefile, modify the compiler definitions as necessary.
d. make # produces hydrodynamics
Please send suggestions, comments, and questions.
Thanks,
Jeffrey D. Oldham
oldham at codesourcery.com
-------------- next part --------------
// Oldham, Jeffrey D.
// 2001Mar27
// Pooma
// Hydrodynamics Kernel Written Using POOMA
#include <iostream>
#include <stdlib.h>
#include <cmath>
#include "Pooma/NewFields.h"
#include "ComputeVolumes.h"
#include "Product.h"
#include "ConstantVectorComponentBC.h"
// This hydrodynamics program implements "The Construction of
// Compatible Hydrodynamics Algorithms Utilizing Conservation of Total
// Energy," by E.J. Caramana, D.E. Burton, M.J. Shashkov, and
// P.P. Whalen, \emph{Journal of Computational Physics}, vol. 146
// (1998), pp. 227-262.
// The code uses fields with two different granularities: coarse and
// fine. Both grid conceptually cover the same domain, but the
// fine-grain one has twice the number of values for each dimension.
// Currently, Pooma does not nicely support operations between fields
// with different granularities. The code occurring between
// `CORNER_FIELD_TEMPORARY' preprocessor values substitutes for this
// lack of support. See the separately distributed proposal for
// details.
// Preprocessor symbols:
// CORNER_FIELD_TEMPORARY: Define this symbol. Include code to deal
// with different granularity fields.
// PSEUDOCODE: Do not define this symbol. Surrounds desired code to
// deal with different granularity fields.
// DEBUG: If defined, print some information about internal program
// values.
// DEBUG2: If defined, print much more information about internal
// program values.
// Sample compilation command:
// g++ -g -DCORNER_FIELD_TEMPORARY -I/nfs/oz/home/oldham/pooma/r2/src/ -I/nfs/oz/home/oldham/pooma/r2/lib/LINUXgcc/ -L. hydrodynamics.cc -lpooma-gcc -o hydrodynamics
#ifdef CORNER_FIELD_TEMPORARY
/** FORWARD DECLARATIONS **/
template<class InputGeometry, class OutputGeometry, class T, class EngineTag>
inline
void
sumAroundCell (const Field<InputGeometry, T, EngineTag> & input,
Field<OutputGeometry, T, EngineTag> & output);
template<class InputGeometry, class OutputGeometry, class T, class EngineTag>
inline
void
sumAroundVertex(const Field<InputGeometry, T, EngineTag> & input,
Field<OutputGeometry, T, EngineTag> & output);
template<class InputGeometry1,
class InputT1, class InputT2,
class InputEngineTag1, class InputEngineTag2,
class OutputGeometry, class OutputT, class OutputEngineTag>
inline
void
computeCornerMasses(const Field<InputGeometry1, InputT1, InputEngineTag1> & pressureDensity,
const Field<OutputGeometry, InputT2, InputEngineTag2> & fineCoordinates,
Field<OutputGeometry, OutputT, OutputEngineTag> & output);
template <class InputGeometry1,
class InputT1, class InputT2,
class InputEngineTag1, class InputEngineTag2,
class OutputGeometry, class OutputT, class OutputEngineTag>
inline
void
computeCornerForces(const Field<InputGeometry1, InputT1, InputEngineTag1> & pressure,
const Field<OutputGeometry, InputT2, InputEngineTag2> & coordinates,
Field<OutputGeometry, OutputT, OutputEngineTag> & output);
template <class InputGeometry1,
class InputT1,
class InputEngineTag1,
class OutputGeometry, class OutputT, class OutputEngineTag>
inline
void
copyVelocities(const Field<InputGeometry1, InputT1, InputEngineTag1> & input1,
const Field<InputGeometry1, InputT1, InputEngineTag1> & input2,
Field<OutputGeometry, OutputT, OutputEngineTag> & output);
#endif // CORNER_FIELD_TEMPORARY
/** THE HYDRODYNAMICS PROGRAM **/
int main(int argc, char *argv[])
{
// Set up the Pooma library.
Pooma::initialize(argc,argv);
#ifdef DEBUG
std::cout << "Hello, world." << std::endl;
#endif // DEBUG
/* DECLARATIONS */
// Values
const double gamma = 4.0/3; // gas pressure constant $\gamma$
const double timeStep = 0.01; // $\Delta t$
unsigned nuIterations = 1; // number of iterations
// Create a simple layout.
const unsigned Dim = 2; // Work in a 2D world.
const unsigned nHorizVerts = 11; // number of horizontal vertices
const unsigned nAngles = 5; // number of angles
Interval<Dim> vertexDomain;
vertexDomain[0] = Interval<1>(nHorizVerts);
vertexDomain[1] = Interval<1>(nAngles);
DomainLayout<Dim> layout(vertexDomain, GuardLayers<2>(1));
// Preparation for Field creation.
Vector<Dim> origin(0.0);
Vector<Dim> spacings(1.0,1.0); // TODO: What does this do?
typedef UniformRectilinear<Dim, double, Cartesian<Dim> > Geometry_t;
typedef Field<Geometry_t, double, Brick> Fields_t;
typedef Field<Geometry_t, double, Brick> ConstFields_t; // TODO: Change to ConstField when ConstField is available.
typedef Field<Geometry_t, Vector<Dim>, Brick> Fieldv_t;
typedef Field<Geometry_t, Vector<Dim>, Brick> ConstFieldv_t; // TODO: Change to ConstField when ConstField is available.
// Cell-centered Fields.
Cell cell;
Fields_t internalEnergy (cell, layout, origin, spacings);
ConstFields_t zoneMass (cell, layout, origin, spacings);
Fields_t pressure (cell, layout, origin, spacings);
Fields_t pressureDensity (cell, layout, origin, spacings);
Fields_t volume (cell, layout, origin, spacings);
// Vertex-centered Fields.
Vert vert;
ConstFields_t pointMass (vert, layout, origin, spacings);
Fieldv_t coordinates (vert, layout, origin, spacings);
Fieldv_t velocity (vert, layout, origin, spacings);
Fieldv_t velocityChange (vert, layout, origin, spacings);
// Corner Fields.
// TODO: How should I implement corner Fields?
#ifdef CORNER_FIELD_TEMPORARY
// For now, implement corner Fields using a Field with twice as many
// entries in each dimension.
Interval<Dim> cornerVertexDomain;
for (unsigned d = 0; d < Dim; ++d)
cornerVertexDomain[d] = Interval<1>(2*vertexDomain[d].min(),2*vertexDomain[d].max());
DomainLayout<Dim> cornerLayout(cornerVertexDomain, GuardLayers<2>(1));
Fieldv_t cornerForces (cell, cornerLayout, origin, spacings);
Fields_t cornerMasses (cell, cornerLayout, origin, spacings);
Fieldv_t cornerCoordinates(vert, cornerLayout, origin, spacings);
Fieldv_t cornerVelocities(vert, cornerLayout, origin, spacings);
Fields_t tmp1 (cell, cornerLayout, origin, spacings);
Fields_t tmp2 (cell, cornerLayout, origin, spacings);
#endif // CORNER_FIELD_TEMPORARY
/* INITIALIZATION */
Iota<Dim>::Iota_t range(coordinates.domain());
#define X (range.comp(0))
#define Y (range.comp(1))
coordinates.comp(0) = 1.0*X/(nHorizVerts-1) * cos(Y*M_PI/(2*(nAngles-1)));
coordinates.comp(1) = 1.0*X/(nHorizVerts-1) * sin(Y*M_PI/(2*(nAngles-1)));
Iota<Dim>::Iota_t crange(cornerCoordinates.domain());
#define X (crange.comp(0))
#define Y (crange.comp(1))
cornerCoordinates.comp(0) =
0.5*X/(nHorizVerts-1) * cos(0.5*Y*M_PI/(2*(nAngles-1)));
cornerCoordinates.comp(1) =
0.5*X/(nHorizVerts-1) * sin(0.5*Y*M_PI/(2*(nAngles-1)));
#ifdef DEBUG
std::cout << "initial coordinates:\n" << coordinates << std::endl;
#endif // DEBUG
internalEnergy = 1.0;
pressureDensity = 1.0;
#ifdef PSEUDOCODE
// This is the desired code:
cornerMasses = replicate<2>(pressureDensity) *
computeVolumes(interpolate<2>(coordinates));
pointMass = total<2>(cornerMasses);
zoneMass = total<2>(cornerMasses);
// This is the code that Pooma currently supports:
#endif // PSEUDOCODE
#ifdef CORNER_FIELD_TEMPORARY
computeCornerMasses(pressureDensity, cornerCoordinates, cornerMasses);
sumAroundVertex(cornerMasses, pointMass);
sumAroundCell(cornerMasses, zoneMass);
#endif // CORNER_FIELD_TEMPORARY
PInsist(min(pointMass) > 0.0, "Zero masses are not supported.");
#ifdef DEBUG2
std::cout << "pointMass:\n" << pointMass << std::endl;
std::cout << "zoneMass:\n" << zoneMass << std::endl;
#endif // DEBUG2
velocity = Vector<Dim>(0.0);
velocityChange.addUpdaters(AllConstantFaceBC<Vector<Dim> >(Vector<Dim>(0.0), false));
for (int d = 0; d < Dim; ++d)
velocityChange.addUpdater(ConstantVectorComponentBC<Vector<Dim>::Element_t>(2*d, d, 0.0, true));
velocityChange.addUpdater(ConstantVectorComponentBC<Vector<Dim>::Element_t>(3, 0, 0.0, true));
/* ITERATIONS */
for (; nuIterations > 0; --nuIterations) {
#ifdef DEBUG
std::cout << "Begin an iteration." << std::endl;
#endif // DEBUG
pressure = (gamma - 1.0) * pressureDensity * internalEnergy;
#ifdef DEBUG2
std::cout << "pressure:\n" << pressure << std::endl;
#endif // DEBUG2
#ifdef PSEUDOCODE
// This is the desired code.
forces = replicate<2>(pressure) * addNormals(coordinates);
velocityChange = (timeStep / pointMass) * total<2>(cornerForces);
// This is the code that Pooma currently supports:
#endif // PSEUDOCODE
#ifdef CORNER_FIELD_TEMPORARY
computeCornerForces(pressure, coordinates, cornerForces);
sumAroundVertex(cornerForces, velocityChange);
velocityChange *= timeStep / pointMass;
#endif // CORNER_FIELD_TEMPORARY
#ifdef DEBUG2
std::cout << "corner forces:\n" << cornerForces << std::endl;
std::cout << "velocity changes:\n" << velocityChange << std::endl;
#endif // DEBUG2
velocityChange.update();
#ifdef PSEUDOCODE
// This is the desired code.
internalEnergy +=
-(timeStep / pointMass) *
total<2>(dot(cornerForces, replicate<2>(velocity + 0.5*velocityChange)));
// This is the code that Pooma currently supports:
#endif // PSEUDOCODE
#ifdef CORNER_FIELD_TEMPORARY
copyVelocities(velocity, velocityChange, cornerVelocities);
tmp1 = dot(cornerForces, cornerVelocities);
sumAroundCell(tmp1, tmp2);
internalEnergy += -(timeStep / pointMass) * tmp2;
#endif // CORNER_FIELD_TEMPORARY
coordinates += (velocity + 0.5*velocityChange) * timeStep;
velocity += velocityChange;
volume = computeVolumes(coordinates);
pressureDensity = zoneMass / volume;
#ifdef DEBUG2
std::cout << "pressure density:\n" << pressureDensity << std::endl;
#endif // DEBUG2
}
/* TERMINATION */
std::cout << "final coordinates:\n" << coordinates << std::endl;
#ifdef DEBUG
std::cout << "Goodbye, world." << std::endl;
#endif // DEBUG
Pooma::finalize();
return EXIT_SUCCESS;
}
/** HELPER FUNCTIONS **/
#ifdef CORNER_FIELD_TEMPORARY
// Following is temporary code to permit compilation using the current
// Pooma code. Adding new Pooma abstractions may eliminate the need
// for this code.
static Loc<2> offset[] = { Loc<2>(0,0), Loc<2>(0,1),
Loc<2>(1,0), Loc<2>(1,1) };
// Given a Field with twice the refinement, form a cell-centered Field
// by summing the 1<<Dim corresponding corners.
template<class InputGeometry, class OutputGeometry, class T, class EngineTag>
inline
void
sumAroundCell(const Field<InputGeometry, T, EngineTag> & input,
Field<OutputGeometry, T, EngineTag> & output)
{
const int dim = Field<OutputGeometry, T, EngineTag>::dimensions;
Field<OutputGeometry,T,EngineTag>::Domain_t range = output.physicalDomain();
output(range) = 0;
for (unsigned counter = 0; counter < (1<<dim); ++counter)
output(range) += input(2*range+offset[counter]);
return;
}
// Given a Field with twice the refinement, form a vertex-centered Field
// by summing the 1<<Dim corresponding corners.
template<class InputGeometry, class OutputGeometry, class T, class EngineTag>
inline
void
sumAroundVertex(const Field<InputGeometry, T, EngineTag> & input,
Field<OutputGeometry, T, EngineTag> & output)
{
const int dim = Field<OutputGeometry, T, EngineTag>::dimensions;
Field<OutputGeometry,T,EngineTag>::Domain_t range = output.physicalDomain();
output(range) = 0;
for (unsigned counter = 0; counter < (1<<dim); ++counter)
output(range) += input(2*range-offset[counter]);
return;
}
// Given a quadrilateral's two diagonals as vectors, return the
// quadrilateral's volume.
template <int Dim, class T, class E>
inline
double
quadrilateralVolume (const Vector<Dim,T,E> &v, const Vector<Dim,T,E> &w)
{
return 0.5 * std::abs (sum (product (v, w)));
}
// Given a pressure density Field and a fine-mesh coordinate Field, compute
// corner masses, placing in "output".
// TODO: Extend to multiple dimensions.
template<class InputGeometry1,
class InputT1, class InputT2,
class InputEngineTag1, class InputEngineTag2,
class OutputGeometry, class OutputT, class OutputEngineTag>
inline
void
computeCornerMasses(const Field<InputGeometry1, InputT1, InputEngineTag1> & pressureDensity,
const Field<OutputGeometry, InputT2, InputEngineTag2> & fineCoordinates,
Field<OutputGeometry, OutputT, OutputEngineTag> & output)
{
// Assume the temporary finer coordinates Field fineCoordinates
// already has values.
// Compute the volumes in fineCoordinates.
output = computeVolumes(fineCoordinates);
// Multiply each volume by the pressure density.
Field<InputGeometry1,InputT1,InputEngineTag1>::Domain_t range = pressureDensity.physicalDomain();
const unsigned corners =
(1 << (Field<InputGeometry1, InputT1, InputEngineTag1>::dimensions));;
for (unsigned c = 0; c < corners; ++c)
output(2*range+offset[c]) *= pressureDensity(range);
return;
}
// Given a pressure Field and a coordinate Field, compute the corner
// forces in a fine-mesh Field.
template <class InputGeometry1,
class InputT1, class InputT2,
class InputEngineTag1, class InputEngineTag2,
class OutputGeometry, class OutputT, class OutputEngineTag>
inline
void
computeCornerForces(const Field<InputGeometry1, InputT1, InputEngineTag1> & pressure,
const Field<OutputGeometry, InputT2, InputEngineTag2> & coordinates,
Field<OutputGeometry, OutputT, OutputEngineTag> & output)
{
Field<InputGeometry1,InputT1,InputEngineTag1>::Domain_t range = pressure.physicalDomain();
output(2*range+offset[0]) =
pressure(range) *
0.5 * product(InputT2(1.0),
coordinates(range+offset[2])
- coordinates(range+offset[1]));
output(2*range+offset[1]) =
pressure(range) *
0.5 * product(InputT2(1.0),
coordinates(range+offset[0])
- coordinates(range+offset[3]));
output(2*range+offset[2]) =
pressure(range) *
0.5 * product(InputT2(1.0),
coordinates(range+offset[3])
- coordinates(range+offset[0]));
output(2*range+offset[3]) =
pressure(range) *
0.5 * product(InputT2(1.0),
coordinates(range+offset[1])
- coordinates(range+offset[2]));
return;
}
template <class InputGeometry, class InputT, class InputEngineTag,
class OutputGeometry, class OutputT, class OutputEngineTag>
inline
void
copyVelocities(const Field<InputGeometry, InputT, InputEngineTag> & input1,
const Field<InputGeometry, InputT, InputEngineTag> & input2,
Field<OutputGeometry, OutputT, OutputEngineTag> & output)
{
Field<InputGeometry,InputT,InputEngineTag>::Domain_t range =
input1.physicalDomain();
output(2*range-offset[0]) =
output(2*range-offset[1]) =
output(2*range-offset[2]) =
output(2*range-offset[3]) = input1(range) + 0.5*input2(range);
return;
}
#endif // CORNER_FIELD_TEMPORARY
-------------- next part --------------
begin 660 hydrodynamics.tgz
M'XL(`(LSUCH``^P]:W<:1[+Y*GY%'R='`7D8'I:E#9*U5P)D*RL+'0G9-R?9
MXQV&!F8]S'!G!A'BX_WMMQ[=\^(A6P8GN3?L)H+I[NKJJNKJ>G1-FOYX,HWD
M&]^=CF5HCK[9PJ=:JU8/#JK?5*O5VN'S??Q;/3PXI+_PJ%X_J']3/:P?[!\<
M'NYCO]JSZL'^-Z*Z#63RGVD868$0W_AN?V2-5_>S0_=KH/.U/Y5*>9.?0J4B
MFJX5AC)LX'<AFAD!PV<O7;]GN>)\ZMF1XWNB*\<3UXKB$79^Q&8Q+'SK#+R^
M'(CK3N?UZ;M7/[5N.JV?KDY?7S1OW]UVVU?-B\MWKPK?0A?'DP_T`NPV\2EL
M?)5`R<Z]#.X=.6L(^)5C#/W,T%F(AA#\O<P<L(@[_D0&5N0'HCA`AOE!2<!3
M`!W)7W&:G1U!'UMZD0QD7YP[TNV'IIJDO,'_(;R'!0C_R8I0L01K>QM8$UB*
M&.A15N!/O7YJ!83X;20]VW&/L\0Y:302.I1,(;HCF1J9HZ0B5"@L.YI:KCL7
MTQ`(8P7I,1,KB!Q813B1-OQU?B-R$]U$ZK-Y\=^"I'7G$PG;A>B_6>#;P/;"
ML]UI7VX#VV\=ABV>7,D9R5.%_FV.GJQLN_"<B`1`!N&*?BUG,.A,PDI:0K-=
MKP._/[4C?+ at -DIW[P<P*^J(E;=<*6%+_%,PFY*>N"YB'=N!,$/-82ZS5ARNT
MI`.[6F]P82,`$:*:]8;""D4T`EZ<<^L3$2F]A&!@NUMC&:$"@H$95:/Z&^U?
M)X$,0\#PQ!"SD6./A`,`Y!@T*\R*4/H.K`*@B'M")U;.I#6NEVN4!(C6VK&J
M=CQ8FK!("X("=]VDR9]&T/:'T^(OG7L)JGMA)4100P!MIH%'M$ITLQ=&U"QF
M3C3*+3.T+1!H`3\#AVP-/4S^#VIN9"OR%`>%WRNJZ[-ML[)*9-XD0"2Q]$`^
M7-'U(_AW6APW.E-!0S[F#=&M&VIKO)0^"'TP/RGP[^QF.BIL9=6K]H$_6&*,
M;G3R;:PFIW_>2-05E3<@_Z)\(FY)?BM-D,\_$EN7C'1`_[2<L0;1?:V_-6^7
M at SGFI1[3H&X==.*=YX#R'-_`<P>VIK1T(\`"*.*D\*$PF?9<QVX4"D)$;)7`
M6-C/'=)G;5:$[Z(C[(!44XJN20H!MCMH&1LU!C1_**`B4"H%^Q9+1_#H(PZ]
M]&V<^X15:&<P"&44+HP%CK"*M35X/#^0=XTT;`VL6$LF0!5(_R;]%4S)H#15
M`_^!M5E3E[0Z'AH(.@"UY:!>Z\WCYI"/&=!Q0_1G[`2>243(*UGQ at 3'(\?"\
M?;+8FX!!D]@=H(9MTQ0E#0&T)"'Z6D8COZ^PQP:4!=>?R:#]*R`7%?%W7Q$/
M!BNR5(\`#'>>HNF^MG.-.Y/)O(C["<%!F<D+`C3H0[18TNN!Y;#`.K45XC`H
M.K6GT%CF;]6$<9N;6OVMKT"A:CX7>\#=?J-A]4)1#*=C49RP'5C8 at 4^1D32<
M.F!:9CSQ1[5D4+M('NGV&K?#9XOK47^?+=DM=Z$4<!*#O:)/9K4 at 5-ZZ4XRI
M`3`,^DEKR/R$'[`J,T4P)1V9T<5J2>RI!E',0(+-N)>!541V%W:RG>J+G4H*
M7GGIA+75$U;SL.J?,F$UGO#IT at GKGS'ALA4ND*&JA>/CEL[O!XW#K9]T^GR+
M#RO]A15<UQJ>%/!T\<"D+Z1M^5NRM7/!@V,8'I^81D$%&HZ32:`]`2Q.&@WT
MIV$_V4N5[=K!NX-2X4/J[%N#"5N%>!0FW9>L1?4R'L):W,B(@2G]1+\;C;'U
M7A85D.+`3)T3Q1)N5Q`FT#/?2J_O#$1E;WV\3>Q5"K]WQ/3_UH?,"\N+V-A"
M<?$]T.9GS0VF`M;'_^OU_><UBO_7#VO[AX<8_W]>>_97_/^K?%"5[I5%\^E3
M_(N:];1YV0`36GKA!CU%BA*.P$(-_4$TLP*)40?X!_0,G/)SL`C0LN>(<W$$
MKKGC"3``7#!EG]QVSKMO3V_:3TH4%."GI"9$$7P]_.V*3N_?(,'E#KCQ8&KW
MM<U)TYQ.)K`@]@)+8"8CF+'5!QSN+<>U>F!M3+V^#,C5!SM]'(+QC.&B'@!"
M9#C@*S0B8F2%HB>EAX"L"9 at H]]`1XSJ!=*4%- at Q%&N#`\,$!16PN3\O-IKB:
MCGLP"_TH__"W\L%S4Y]W=YXKX7SQ`8- at Y@`$!_2AC6,-P`KHMC"WL*:P0F7K
M6X0*'F+^7$H!F.COY/;BNL!Y(D<AFI,C#'XQ(.PYEI%$^+DG0KKT0W'J6F/X
M<T5T@[/XTNKYU&^NJ`6J(PHLL,VN?%.\+1_N5Y^7VU<OR\\.:/T(AR8V;TW1
MDACOIB`4S-+V9#"<*ZI2^TL at 8>!1.RXP<(:C*!21CY%S`R$%DBU!:1!+^PXX
M,4X/CK8L>4P"R4X@\'@.)N8$SJJD.P&;!(B.!"X'SCVL[QXXY@?O65IXL#O'
M06`4S'/TQZ6!S\@!JI$5#`$CE`"G#ZR(1E;$_:_\"/:0DO(YF.G`2[U\YEPX
M<B:8$Q"9Y?4QNP(BC8@[, at 0:74D'I8)HF2*3YP=YOK)<OY<A30F[+`#%#JN7
M'%A$L<`SW4&A at N\@GQ12P,ZN8_7`IT;A"!1"(9P$H:.>#M1LP`X6DS3-M1!?
M#!(R0?O8AR/=0:KXRHQ?I+ at APJD]BKLBE!A$"'1V0;S!9X5=%0!'QE;P'I$/
M?0RV>GZ$H,$L&B!:3L3;3HFP]HJ3/3X(_#%LOJO+&.%S6-48]E!&_8"8`WM)
MO\!,8)<(61X##%J&[X^M_[)LUW0MSS6'_CW)$X"Y=X`%1"+63#/9$Q-K"!(`
M$AU%DT:E,IO-S/30"D&KF)O5LEGUO1FX*?`4,X]3NVP[G%NV/&MN.U%[U7Y[
M?M&^;+V[NVZ==MLW[YJ=J]ONZ57W3;O9[=PT.Z^O.U?MJ^Y9<R%[^UE#MYZ*
M%2MMKH8X@[9^$W3%$/6LSB]$D0I$ZP`UA39I;]U+SD1H("#*I#Z&4\S7@/J"
M;0!;:.ZC(LKE%NF3@"JZ_M#!XW5>@NT?Q]@$Z(B1&`"/<0=8[+0#]F(9M'O+
MG4KSSY`9^BMIJ8"GLHAWR/,(CIW*E8]IR6R.,6F]/@U#&7"V<4FZ\F[2!PD.
M]-\52<U<KS.PF_[DZ<L5403PC]<$$M0#H at D'5#:+U=H$Z"I%Q&G.;`(%SWW%
M*ZV at DI,_(@M&AIRJ%#V\6H%V/2B+OL.C,6T7:YBT at J'$)IR$H+EF:&W3*4J*
M!*<$C!:5'``C4XNTDJ.2HUHQL0T(6E,OKUA*;GTH4PV1'$O+(QMS9($Y0B`B
M9RS+?3F!4Q^G`2]8B--TR%R$8&.`P0)K\&><[P5J`0R)%H0MV>BUV!)`Q5F.
M2<$K`F1_DX%_Q":4-1C`TL`R\]RY,EH6TYH\,JS$RT>%;EMH\^`L\03LV/#!
M1T97CSN`?.G4;\PPI%`WCNVFLM,.+PH#0MISB!7^MI1[C(9.%R4)H17RF4WT
M;!2A)5D7C#P7Q&IL*"6!HLB1[?ML.P7;Z884Q>^Z8E<3-&[I^;ZK)4C/<J;9
M^@)`NZ$L(0H-FN;=N(A_2@;%=7/306/NB>JG9X4.^FO)X*#ZBJFAYXJ6DN"!
M'U3292UU>-VK=<TN6/_2S:Z/'IG\0V.YN%#NM?!\<<'<,7GP\,)YQ,KV4KS\
M+0D at G+#.D#V]^)*':/D20Q-R at KZ;:C8?$$ZQJ\>_>)@3P2 at LQ?D8IC[('SQ5
MK#A:S at C59^'Y48X/JE_RX&@]%U3_E>V4Q8T#SGNH4E7>:DM,L6WPHT$E5)3#
M"K8U>L9P*.)Q at M&=4/%#JX3B0J9243+.;>:HMCA at D:PTMAL3=G%(BL+4=YV"
M61R]FMYQ;OH.[#^*K\SP8$`O?&1Y0YDY+[X/^<0PA#2'IL$QA/0!:TO*0-_[
M3C]S7'>%)V?Z%V:3,Q*4:N/%X?\G%%:0C4**]"Q>2TC,#5V1E<-U-,(>VTIW
MG7KJGN at RJXG/0:3<ZK/P:V7#UMNQZLA6%N)Q-EM$0S,)(T/G)-5GY?*.=:)-
M=!N-.->,"2>\0Z$B?2D?XK-G?C0"RV^;/)8`CR0"0<'(>C:1MP&"/!H;G'1)
M6G$-!N)V2N%['L7_W[91ASQ35%+G8HR$V`W!--,&VV>38;=GE\B684(40T/`
M$Z8R^SE at 7US[\*W18#^[J'+X'Q0G=+($W0(5T0TQH"<<MLA#QI1"Y9Z//\L#
M=>M>\4[]58YZ40THEDQO.H9U\AU]T/PO7 at B^NY(:A#%8B7=N:"[&V,QVT<L`
M=9R`CC#)VZ*6XA*8+R5[17V'?#_?,PNJD6[T`*B>;:HCLR+J1X4T,22G1I+!
M0)&)Y:6I\3WZ<GC/46$L,N,#"4H4CJ*9%/^>`E?MD3_!NU\$(1TL,T1DO<=P
MFXVQ>3C9$C!PQ/5==<91EH>\58"(/?GDH'O`Y.$EVCP!4*1D$9Q5Z.P0<U&S
MHVM6CF^(E7*D1N)8?<19DQ2$H`BD6GV<[^Z*>$LG[,EDWDU?7UG[N?]/E():
M*1G#TP%#JGI*"0[(8G--[=98:%P42B3/R!F.,"H.3C(=Q0FK&7?%YMW,M!_B
M;P3N%<+ at P2+;HB7)XWP:N*F9:*?#8D%WQI9(FR:J]Q('A1D1SM*(H%X24/A%
M\(K]6+)3/%)XG3)A%G=-YI-ZIK<1\4!<H`2`Q71<.RFF6LRQ]6L1;YIIA)\J
M#H`/LZ.O*R]^%B"D$/\H4GQ=R8#+F'G+E_MI;*![?AMD`\';`!MV/I\/J->8
MZ#L+%,8VX(Q>35GQ*$WT`O][2[Y)2V:N=?Y'GVVE;7JI\<U.Y1]-)NX\&Y2*
M35GM&)&U/R7DBME3+^9_3/62B3&OHJW"C<#]!6<)CK#"3JI#XA,=%1*J%S(.
M0GQNQR=^H\'G%IS]>N[8X-\@T93![[H/!%ZS]2443T1E38D=S!^S004'WEDS
M5(&^M>'<+<;L\B&[=8O+&LSK>A:7Q='6!-#$)\;0R"S[A(`93Y6)E'&DPNZM
M"XI]T/[HPQY4RG%2>R)U=_9!>WEWD%R>33;0X\SME`)%:X1([]#A#W^.15WL
MI?=)'QQX#Q/K(30_?5J*U=O`M/I]K716V\S=DZ)C++*`;R3O:*,;'C"M2_$V
M7G#QUSGU:8\>P?!FUI<,*X])"__>%[6V](D+![<XQP/U_]7]6I7N_SW??U:M
MU0^@?VW_L/;7_;^O\=G"^?90?;^N4?BJA?W7\/VNV7V at L#_I5?@S%?8_7$JI
M:([E/5CG/[;`2,/[3C8,4]E5/G]@>.Y>QY_B>L1?ESD8>.JB1=?QYA4^@+=T
MKV+MW8+/D<F&N)8!WL##K/PZT32%N``+G:[^O??\&8)QP#!&:]X0O6F$MP#!
MA)Z&<C!U,4#CJ2IB7<FK"Z<'&+2QPJU<5>+ at W9L<[G3U0&K<&-O,8MB9M@@[
M?7G6]@.^B=G'!61HDQ"%INRT. at W15%D8/W]W`RV[EC-^43=3V?9\269LFYX4
M4K6770-^:UZQ"9YK%;OWVCA?:)EQ/8IZ7N>GEA>"$X\F&G_#,JP7XIZJL<1,
ME4&JIAHVE>^I<`K:.%JIZ[P5G(];<=26BG=*:M4-4^6@,3LPK"?^I=PS?=MB
MLX at 5V,\7YYZRW8C`U^UN^UW[]77W)[9B;^[0A+TMQKU*BVY)MQ:SOIXJZXO=
MXS/'`[?FAFA]C)VQCC<&F"I46E+XUZV)72MVV>IBM[>BC%&7+!8MHY>N+,N[
M4#4C1EI):<T@`:[K!OVE#8O)+^6U]5ZJA?"JKOR^/([78J2C`DU0#)&\E-9`
MITYJ!LS:KE%M%CY_MW+`:1!8\V/`";`!/)(!]#5YPP,0+;NE\C/M"M?@ECS(
M71%D:[Q6+>@12_KL174#*3-57RE"<]N)+@&C7T6:9#TZU-U5'O=Z-*AK0.6J
M'S]-9!)Y_[U$!#;#!J0"M]17%X0LZE^3]SSS"G:3&EC-[)?;U!!**#^=_6K1
MGZ\A\C,ELI`'^86"L79)G[VH+Y:2)>BLDI(E:*S5$-VO(B)P%&Y/*O"<W98@
MY!'_^KQG#![-[O55Q;&S__^YJG at T[X,8S(&ICAV:MKV-.=;'_YX]K^W7=/RO
M5GN.\;]G,.*O^-_7^&`XB59NB!_E8!#(N6B1WU('9KVV at OHAO5+)]\<6N;:O
MTO(B_B$#3[KB;>!$$;C:=R%ZJ;374I&(8\<'ET5:XY/4LS#JNT[/'*6?V>C?
MGJ3?:H?35G3%2.[M>-E7&ZQ\'5YZP,I2]R>)VY[9$.CM#0-KG'X1VQ.\C1)?
M)>+""/VB)G"VL=(P2Z13=^@#@4;C4%`)S6]()`2`669=64%O2$`P7!MK/,&R
MWK;YHRF:>#??\L"9:IEM4YQ-@\CW#/$:VVY'5CAZ[]\;^B5NU^:U*=Z.+%="
MEU_ at M!M]^-$'=6NYR>NO=$'O]6@>`GX?#0R*F**V?X`0BK4??OA;R1"3B2GJ
M]<-R_:!N%G3)@`TG!P8P0L'WC53%Y<P7H&<',L#KRD`O;^I:`94*-6"(%812
MXX>15U.(,Q^�.GCSZA+2?J-9DVEK?RG1X\9#@5RY$25=6)X\LP@^-180A>
M@(IF6&B;O8C`%1,4]9"6/1)QUHI$NSD-$%-W;K!@PT22RTD]``5XA-/)Q`^B
M),2"%1?1#*NN>=T(A9:^8MFJM)G(Y=LVS(<\5S!P\+^:G9NK]LT[SC]UP7'O
MW)S>_/0]E28'/E\SULL(I[TP<B(0]E`5XG*IAVO9[W&U"EV8]%8R)4*L;X9C
MDJN85>ES'Z5YXH?`>H#"M2*1Y;@A\_<Z/7,X'_=\EV-URU%MB!;'T:F0A?MC
M>$[M-EIZ!'(A6:PY4;N29G.A+[`A)K?MNU:GV6FU81*.C_67S74+=,7<*Y7K
MTVNOU*Q8)JFO=L#\GSAMJWUV][*!5<RJA at 8+N]&W"/WQLA)A!^^)X$922H)G
M99XE`.O+((ZYZ'EIY;$&JQ:A-9"&2R72%JHCBFPZ+H^%[Z`D^L2P(;[,82C*
MK>6<$^6+BC<(*_YOE1&LK,(GGZI&#NJ5,+`K#_4!W5VYO+BZ^^^AC9TO39$W
M)$39I=[E(7[WL^V4JZ'WTRS%$!:YMR?..S=O3V]:0,;FY>G-:?>B<W4K]L!6
MRYOQ%_B.M?P[>_CE5)_R)A\V\PN8A"^$T_$IR12]&R[CA>=FR27D^45O=&V4
M/CPFCT1^$+]\KK08[MK&DMY0"=ACUK2SR=744O>-4\U),))_UI=WBR?,=H\?
M+PY;3C5^VLW^7$%!]1ZFI at _&3O":WAF[FHJ`5[R>',9((W+DIH%LX7D4S9&X
MA&H*7!YA38_\2A$<:I6F[P=]<.S@>-#@E at .*UYQ?[5+>B8>8]ZG<>RS[ML._
M<[RZ\^7\^U+&V7\,ICW(I:VP8S)_(UW?)E/IT9P at Y51+-.X7P:D_I+D_B0W)
M!:-5YQJ?;-U7;9&)18CKF\[+F]/7?+R1A8"WZ.GN=S"T#7K3C-B#[_<__Y-#
M.ME;^FS'PJD<6%B6+82J+W#BEXL7"1!"0$3Y]"4#!3K32R1M-#Z.C\635W#T
M^0:^G\7MFT_P$;7#VMST$GDP8K*7.Z(KZCKH&S)9"IHU?7^*SM'0`DM%O!#[
M9K7R[&AG!PT6*XPW5E)N_-TOU/6[/`"L(;N-Y`3OJ)G5&H/X[I>6="-+1-A]
MZF%U)-B#WO0BBDUXO+!.?1,_P8E;%<H<7<+:=_+X\"ZS/Z4:2T8AAMQRL/BA
MS@#?^L%[>O&&J+<4W19'>*_`"_P-#V%"!7#)H#+"5M^+*.L=X,N#PB4P3KVA
M*W'\\_Q2+&HIB.3Z,KV?ENN^^78K!N+2OW^NYJ\[IW`L+?2N+?1F;*@G][DD
M<O',3+IB&H0A7B9WN8_K^,+;TI&B_'7`/HO.5G/INXW\(+]-I'++)P+0A`U8
M!/[3[.FF<&)A&498K)E5HX8=N`:$$N1O,55+[AZY$GW_[ZGHZ*KW"[/<&>B)
M@[IV+(\G.HGO=:ZK=L)7)VH`]#D+'/O]B?K/9#QF9/)&<QR>K"U)_J??>8Y7
M(%*_\7ZQ?B.2N7;J%$V-#-+W#R&];&2"POV7(JUV:N8][OH_.J)>ZHSO:T<<
M-95BIXJC*P(L>^AA*!G]W_:.M3>-'/B9_156DDK+8PE+>JI$0^]20EM.(:U*
M<J6]RT4)T!PZ"#E(6D55__MY'O;:WA=$3=H[L5+396V/Q^/Q>#SVC"O,3A7-
M.\A4%IV%')LCT/W$<H5U.277"G<MR)KB:H7IF,MR=;*P3 at H]L:!Q-Z?X$7&B
M7,TD88DJ.-"704_VOZG]((ZKE/U$JL/MG0LRRRU1F#D-]4:#()IY7\T^JRAE
M'>.:A8%9X.><I2Z"@T!DE[//E30 at XF9!GF7&109D]Y(SIUSWWQ(8OL8`IB+7
M[.7."P3^-V=VP`/@>K89TBGPH=B%Z>ZI*)>'=-X[7CCN.5,O?;+3V7\FX;OV
M3HI/(E03??'CU29/*!:C13J_,1I,L'E#:6`L^\3J("(LC'4:L][J,`R]>340
MT)+KZ558N",59-GZ2F7S=6'4'#N'G:/.WD'G`RJ/K#MV9M=GV/^-!KS*ZN<P
M9'U#;%2'[&B*:*HSO7WA8T[R&Z+0V"KMO95&<;I-@%P$=+-JK=3?-I6A`$^?
M#68+_WVI>_JFLRW96ZD_,HV0B`$+LX#)`9T%+$Z"@:*!PTNIE!ADD6*00`L7
ML*8(#OI:]:=4JD!:)F4208=YH!>HXF6`SES'\-K'[)?&'Y>XH#$^Y:QOA*LZ
M8(]"H]SYF1,8I<B&[47^W!Q<R;16-S1Y6,(T(9(6AJL=@4!S:E'!W)W0X(CC
MU6S"A8SF40=$\W23G*0I5U0MYM(J3E8FIRUD<0>MFI:@`[6WHG8E%@W3[R:B
M2MZLF&1T<XUG<<:JB!C"<?MGE*$24<;.B]=]V#D5=9:1;7+]W9$H+JY]F/2B
M.N2:02Y9*F+C`WB(3ZG'P7L<]CB88".YZ'98NQ[C;0U2<734PS8_NP55*U0Y
MW>>9PZ".:T%2GR1W&`J^7H+9VI7A#+;P#4\_"A6Z:Z[7GODNN`I[ZQ6U2J(B
M!"1J(ZDU9[BA&56:L1NDA!G*15>%>NEZ?L-(?,LJ=B1TIP*>"H_:M at D%F_[4
M-F(\`QH$@?FM*+ZX at M#M]>>C"S!/7$8&CUS#CA"1B).4]\EL$X"0`]'L2K^2
M(R?C#!QC88:@.5A5EL.)20(V4\22^S>Z;*:*5VB2[-5#V(N;+"P!^C2!R8`@
MV at RU;0 at 1"<81G:1[,I2'$IZ.^&2;M[9=6]9GD8!GLLRD+!6'%LD$*C5%$H&6
MD)UY?,,+(^K/:$*/$(Q)0!>"EF04M$H#<9J0*Q%C<D&YM2=K`4LPJ:-KE)ML
M?@_2N(W3-<\-9]=.9UGLKEM>1IW+Z<CB=^!2:RM`X>.R6"6V^F%,844CQV*\
MU2GY[?D=2E=P9</),?*G$QZ++:4)6`L#@)G3"P";Z[1'%I2U\W(ZN30U7:4P
M+L3B*JN>_K<9S/)R6[(O0G'EMP*>.WR^\KPGI[UNY]!<_]D5RB7+REI\UKK@
MY6PV/+\=+;O#H?=2$`_823%=D-K]SM%I[[C5:O=Z[(A4*HE7[8,W[;?BQ?%A
MRSBVD#T>T`P$867YFCC8TIO-*88&G:6Y&LVGXVOKV`?9A7"4TJC4I_5(J at BQ
M-T37L<O19_Y^=K[`:QM0G8`0?J/)>(I4I$-4(PK#KPX:L72""PS&`[RJKOY,
M4/2 at W\'F\X6_^;4*:&[Z1UB,HC/PQ]#,$<H<`AQ]P%>0;_6,6;<`H3FN6:<8
M`8)]!`=QVR^`.0>)-)TJBH2[H'PY+G0D&!;5ASI#<O<C)"N=N<`M0:H*=67<
MG5H"@A70P4NI]*AB%M"14N:L"?%5L5=XG-".`$9)9'<IJH!2MI%1#DOH1U+N
MU8]=X<O>DZ at 50<_GKZ3MVR"E2$2J2<4=OY29,[G(B3%4E9O at RMR6&)3Y1^2W
MASK at LV8V8K9@!6;[Y^9L.`>I#1?R0J at Z/#-[=@$G<1=@Q2=7YX6^PY@/O+KE
M:)9>THN7N80V$3T+%"D)XLX.O4M<0BD5'3G'*E\,DQBN"L%#$2_+H2.^>.=$
M-+^K^YU9PZ%H[[@,F"H[S at 0MS[#UL4&]MF'X1>/5H1 at 3?7HSN<;[)2-6S!Z;
MZ^-J/_!Q-7429@^O,N+[NY3:!)7-+261V,Q!@B"<3>:CL^$MGBC71VUY^P_)
M:40/P`TV!TI5"Z"X(NYDU1N+7>+%6]JL8S7^G`*HN4.DJN6ETWVJ]]S.2Y"=
M3F<F"E'GP`E/7VRG]T/0E?T5V<@4^U(8)$AEEL="K>,7*(@M$>S.[R=%,#$X
M+6)!C?1-D<.:KJ;(214T*@:#Q(EUX@'UOBFFL$AU?6[R/W=N$J7'MQI462J)
MP[RU$[7QI0HK#:.DML/0PDI^GDP-.,Y4I"8;-/`MP/63(F8(TK.$)\4,S,+[
MPJR6C]E.)F;U^\)L)Q^S6B9F._>%69B/69TQ,X5=IARRI4VR+/F>9WXK&C,;
MIY4._.8"^1:G?=.$!\L.1W0D2`ZV.$*;EI`=@2,[DE/#S-1Z9BKR,2.D5SQH
MH"2:Z<G58K=<`^CW=G!=/YD/>.1_'$]&]UD'^G_7TOR_:^%CY?\=AH]K81WB
M/^[4GZS]OQ_BV=S<3'(`A\_L`5Y[@K_(!5R^X:\6VJ!'>&_T&W9_]-ZU7ASL
MO>PU"\$[B%0<O+N<!:#B!Z!(GLU'GDHO!!=BR^?L1?&'G.:"?;2VISLM8J;.
MUI=7K[OMK^R(B']KRF$Q,<UR5/1:+5FYRB>_!&/8O9],X'UG^WQ\"2]>J]]7
M#=D_[K5/#SK/7_;[IYU#":C=DY at SXI`Q%UZY[!WL*W`'U<*FV.OUCKOMN-L&
M&/TY:#H;\SEL.MZT)Y%`RD5>E9YL3N_U\=N6Q*D9\[_TO-?/?VVWCB!-8JQR
M-A[)M.:CZJPH585'#2'?O(),[_>+`O_C/MGR&6OY^B?\DM47P8]SZQ<N)ZM8
MN2!<Y'O9\`KSJ0@^BE)U)O^-Q_+/PO/^DNOOT3RZ'`B\Q[T">H]VNQ'1S<;@
M#_4NYR))V_U]08`$2#5A0H-T]N^`,'J\=[)Q`3%#1/"W2$!@H[J>O];/^ED_
3ZV?]K)_U\W][_ at 7#@K2A`*``````
`
end
-------------- next part --------------
*** Product.h Thu Apr 12 16:07:14 2001
--- Product.h Thu Apr 12 15:46:46 2001
*************** product(const Vector<Dim,T,E> &v, const
*** 50,53 ****
--- 50,127 ----
return answer;
}
+
+ //-----------------------------------------------------------------------------
+ //
+ // Full Description: product() versions that operate on `Field's.
+ //
+ //-----------------------------------------------------------------------------
+
+ struct FnProduct
+ {
+ PETE_EMPTY_CONSTRUCTORS(FnProduct)
+ template<class T1, class T2>
+ inline typename BinaryReturn<T1, T2, FnProduct >::Type_t
+ operator()(const T1 &a, const T2 &b) const
+ {
+ return (product(a,b));
+ }
+ };
+
+ template<class G1,class T1,class E1,int D2,class T2,class E2>
+ inline typename MakeReturn<BinaryNode<FnProduct,
+ typename CreateLeaf<Field<G1,T1,E1> >::Leaf_t,
+ typename CreateLeaf<Array<D2,T2,E2> >::Leaf_t> >::Expression_t
+ product(const Field<G1,T1,E1> & l,const Array<D2,T2,E2> & r)
+ {
+ typedef BinaryNode<FnProduct,
+ typename CreateLeaf<Field<G1,T1,E1> >::Leaf_t,
+ typename CreateLeaf<Array<D2,T2,E2> >::Leaf_t> Tree_t;
+ return MakeReturn<Tree_t>::make(Tree_t(
+ CreateLeaf<Field<G1,T1,E1> >::make(l),
+ CreateLeaf<Array<D2,T2,E2> >::make(r)));
+ }
+
+ template<class G1,class T1,class E1,class T2>
+ inline typename MakeReturn<BinaryNode<FnProduct,
+ typename CreateLeaf<Field<G1,T1,E1> >::Leaf_t,
+ typename CreateLeaf<T2 >::Leaf_t> >::Expression_t
+ product(const Field<G1,T1,E1> & l,const T2 & r)
+ {
+ typedef BinaryNode<FnProduct,
+ typename CreateLeaf<Field<G1,T1,E1> >::Leaf_t,
+ typename CreateLeaf<T2 >::Leaf_t> Tree_t;
+ return MakeReturn<Tree_t>::make(Tree_t(
+ CreateLeaf<Field<G1,T1,E1> >::make(l),
+ CreateLeaf<T2 >::make(r)));
+ }
+
+ template<int D1,class T1,class E1,class G2,class T2,class E2>
+ inline typename MakeReturn<BinaryNode<FnProduct,
+ typename CreateLeaf<Array<D1,T1,E1> >::Leaf_t,
+ typename CreateLeaf<Field<G2,T2,E2> >::Leaf_t> >::Expression_t
+ product(const Array<D1,T1,E1> & l,const Field<G2,T2,E2> & r)
+ {
+ typedef BinaryNode<FnProduct,
+ typename CreateLeaf<Array<D1,T1,E1> >::Leaf_t,
+ typename CreateLeaf<Field<G2,T2,E2> >::Leaf_t> Tree_t;
+ return MakeReturn<Tree_t>::make(Tree_t(
+ CreateLeaf<Array<D1,T1,E1> >::make(l),
+ CreateLeaf<Field<G2,T2,E2> >::make(r)));
+ }
+
+ template<class T1,class G2,class T2,class E2>
+ inline typename MakeReturn<BinaryNode<FnProduct,
+ typename CreateLeaf<T1 >::Leaf_t,
+ typename CreateLeaf<Field<G2,T2,E2> >::Leaf_t> >::Expression_t
+ product(const T1 & l,const Field<G2,T2,E2> & r)
+ {
+ typedef BinaryNode<FnProduct,
+ typename CreateLeaf<T1 >::Leaf_t,
+ typename CreateLeaf<Field<G2,T2,E2> >::Leaf_t> Tree_t;
+ return MakeReturn<Tree_t>::make(Tree_t(
+ CreateLeaf<T1 >::make(l),
+ CreateLeaf<Field<G2,T2,E2> >::make(r)));
+ }
+
#endif /* POOMA_HYDRODYNAMICS_PRODUCT_H */
*** hydrodynamics.cc Thu Apr 12 15:51:54 2001
--- hydrodynamics.cc Thu Apr 12 15:55:42 2001
*************** int main(int argc, char *argv[])
*** 158,174 ****
/* INITIALIZATION */
! // TODO: Is coordinates initialization necessary? What is the best way to do this?
! for (unsigned xIndex = 0; xIndex <= vertexDomain[0].last (); ++xIndex)
! for (unsigned yIndex = 0; yIndex <= vertexDomain[1].last (); ++yIndex)
! coordinates(xIndex,yIndex) =
! Vector<2>(static_cast<double>(xIndex)/(nHorizVerts-1) * cos(yIndex*M_PI/(2*(nAngles-1))),
! static_cast<double>(xIndex)/(nHorizVerts-1) * sin(yIndex*M_PI/(2*(nAngles-1))));
! for (unsigned xIndex = 0; xIndex <= cornerVertexDomain[0].last (); ++xIndex)
! for (unsigned yIndex = 0; yIndex <= cornerVertexDomain[1].last (); ++yIndex)
! cornerCoordinates(xIndex,yIndex) =
! Vector<2>(0.5*xIndex/(nHorizVerts-1) * cos(0.5*yIndex*M_PI/(2*(nAngles-1))),
! 0.5*xIndex/(nHorizVerts-1) * sin(0.5*yIndex*M_PI/(2*(nAngles-1))));
#ifdef DEBUG
std::cout << "initial coordinates:\n" << coordinates << std::endl;
#endif // DEBUG
--- 158,175 ----
/* INITIALIZATION */
! Iota<Dim>::Iota_t range(coordinates.domain());
! #define X (range.comp(0))
! #define Y (range.comp(1))
! coordinates.comp(0) = 1.0*X/(nHorizVerts-1) * cos(Y*M_PI/(2*(nAngles-1)));
! coordinates.comp(1) = 1.0*X/(nHorizVerts-1) * sin(Y*M_PI/(2*(nAngles-1)));
! Iota<Dim>::Iota_t crange(cornerCoordinates.domain());
! #define X (crange.comp(0))
! #define Y (crange.comp(1))
! cornerCoordinates.comp(0) =
! 0.5*X/(nHorizVerts-1) * cos(0.5*Y*M_PI/(2*(nAngles-1)));
! cornerCoordinates.comp(1) =
! 0.5*X/(nHorizVerts-1) * sin(0.5*Y*M_PI/(2*(nAngles-1)));
#ifdef DEBUG
std::cout << "initial coordinates:\n" << coordinates << std::endl;
#endif // DEBUG
*************** computeCornerForces(const Field<InputGeo
*** 354,390 ****
Field<OutputGeometry, OutputT, OutputEngineTag> & output)
{
Field<InputGeometry1,InputT1,InputEngineTag1>::Domain_t range = pressure.physicalDomain();
- for (unsigned xIndex = 0; xIndex <= range.unwrap()[0].last(); ++xIndex)
- for (unsigned yIndex = 0; yIndex <= range.unwrap()[1].last(); ++yIndex) {
- Loc<2> range (xIndex, yIndex);
- output(2*range + offset[0]) =
- pressure(range) *
- 0.5 * product(InputT2(1.0),
- coordinates(range+offset[2])
- - coordinates(range+offset[1]));
- output(2*range+offset[1]) =
- pressure(range) *
- 0.5 * product(InputT2(1.0),
- coordinates(range+offset[0])
- - coordinates(range+offset[3]));
- output(2*range+offset[2]) =
- pressure(range) *
- 0.5 * product(InputT2(1.0),
- coordinates(range+offset[3])
- - coordinates(range+offset[0]));
- output(2*range+offset[3]) =
- pressure(range) *
- 0.5 * product(InputT2(1.0),
- coordinates(range+offset[1])
- - coordinates(range+offset[2]));
- }
- return;
-
- #ifdef PSEUDOCODE
- // TODO: I really wanted to use "range" in a data-parallel
- // statement, but I do not know how to extend "product" to work
- // on Fields for such a statement.
- Field<InputGeometry1,InputT1,InputEngineTag1>::Domain_t range = pressure.physicalDomain();
output(2*range+offset[0]) =
pressure(range) *
0.5 * product(InputT2(1.0),
--- 355,360 ----
*************** computeCornerForces(const Field<InputGeo
*** 405,411 ****
0.5 * product(InputT2(1.0),
coordinates(range+offset[1])
- coordinates(range+offset[2]));
! #endif // PSEUDOCODE
}
--- 375,381 ----
0.5 * product(InputT2(1.0),
coordinates(range+offset[1])
- coordinates(range+offset[2]));
! return;
}
More information about the pooma-dev
mailing list