6 #beginMacro getTravelingWaveSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
12 const int v1c = parameters.dbase.get<
int >(
"v1c");
13 const int v2c = parameters.dbase.get<
int >(
"v2c");
14 const int v3c = parameters.dbase.get<
int >(
"v3c");
17 const int s11c = parameters.dbase.get<
int >(
"s11c");
18 const int s12c = parameters.dbase.get<
int >(
"s12c");
19 const int s13c = parameters.dbase.get<
int >(
"s13c");
20 const int s21c = parameters.dbase.get<
int >(
"s21c");
21 const int s22c = parameters.dbase.get<
int >(
"s22c");
22 const int s23c = parameters.dbase.get<
int >(
"s23c");
23 const int s31c = parameters.dbase.get<
int >(
"s31c");
24 const int s32c = parameters.dbase.get<
int >(
"s32c");
25 const int s33c = parameters.dbase.get<
int >(
"s33c");
27 const int pc = parameters.dbase.get<
int >(
"pc");
37 std::vector<real> &
twd = parameters.dbase.get<std::vector<real> >(
"travelingWaveData");
38 const int np = int(twd[0]);
39 const int ns = int(twd[1]);
43 printF(
"\n\n **************** FIX ME: travelingWave: finish me for HEMP **********\n\n");
50 if(
mg.numberOfDimensions()==2 )
55 real x0 =
X(i1,i2,i3,0);
56 real y0 =
X(i1,i2,i3,1);
60 real s11=0., s12=0., s22=0.;
67 for (
int n=0;
n<
np;
n++ )
69 real ap = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
71 real xi = k1*(x0-xa) + k2*(y0-ya) - cp*t;
80 s12+= -ap*( 2.*
mu*k1*k2 );
85 for (
int n=0;
n<
ns;
n++ )
87 real as = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
89 real xi = k1*(x0-xa) + k2*(y0-ya) - cs*t;
97 s11+= -as*(-2.*
mu*k1*k2 );
98 s12+= -as*(
mu*(k1*k1-k2*k2) );
99 s22+= -as*( 2.*
mu*k1*k2 );
114 state0Local(i1,i2,i3,1) = 1.0;
116 state0Local(i1,i2,i3,0) = state0Local(i1,i2,i3,1)*det(i1,i2,i3)*
mg.gridSpacing(axis1)*
mg.gridSpacing(axis2);
118 i1 > I1Base && i1 < I1Bound &&
119 i2 > I2Base && i2 < I2Bound )
121 real area = 0.25*(det(i1,i2,i3)+det(i1+1,i2,i3)+det(i1+1,i2+1,i3)+det(i1,i2+1,i3));
122 state0Local(i1,i2,i3,0) = area*
mg.gridSpacing(axis1)*
mg.gridSpacing(axis2);
132 if( assignVelocities )
134 U(i1,i2,i3,v1c) = v1;
135 U(i1,i2,i3,v2c) = v2;
139 U(i1,i2,i3,s11c) =s11;
140 U(i1,i2,i3,s12c) =s12;
141 U(i1,i2,i3,s21c) =s12;
142 U(i1,i2,i3,s22c) =s22;
145 real press = -(
lambda+2.0*
mu/3.0)*div;
146 U(i1,i2,i3,pc) = press;
147 U(i1,i2,i3,s11c) += press;
148 U(i1,i2,i3,s22c) += press;
157 ERR(i1,i2,i3,u1c) =
U(i1,i2,i3,u1c) - u1;
158 ERR(i1,i2,i3,u2c) =
U(i1,i2,i3,u2c) - u2;
159 ERR(i1,i2,i3,
uc) =
U(i1,i2,i3,
uc) - x0;
160 ERR(i1,i2,i3,
vc) =
U(i1,i2,i3,
vc) - y0;
164 ERR(i1,i2,i3,
uc) =
U(i1,i2,i3,
uc) - u1;
165 ERR(i1,i2,i3,
vc) =
U(i1,i2,i3,
vc) - u2;
168 if( assignVelocities )
170 ERR(i1,i2,i3,v1c) =
U(i1,i2,i3,v1c) - v1;
171 ERR(i1,i2,i3,v2c) =
U(i1,i2,i3,v2c) - v2;
175 ERR(i1,i2,i3,s11c) =
U(i1,i2,i3,s11c) -s11;
176 ERR(i1,i2,i3,s12c) =
U(i1,i2,i3,s12c) -s12;
177 ERR(i1,i2,i3,s21c) =
U(i1,i2,i3,s21c) -s12;
178 ERR(i1,i2,i3,s22c) =
U(i1,i2,i3,s22c) -s22;
181 real press = -(
lambda+2.0*
mu/3.0)*div;
182 ERR(i1,i2,i3,s11c) -= press;
183 ERR(i1,i2,i3,s22c) -= press;
184 ERR(i1,i2,i3,pc) =
U(i1,i2,i3,pc) - press;
206 #beginMacro getPlaneTravelingWaveSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
212 const int v1c = parameters.dbase.get<
int >(
"v1c");
213 const int v2c = parameters.dbase.get<
int >(
"v2c");
214 const int v3c = parameters.dbase.get<
int >(
"v3c");
216 bool assignVelocities= v1c>=0 ;
217 const int s11c = parameters.dbase.get<
int >(
"s11c");
218 const int s12c = parameters.dbase.get<
int >(
"s12c");
219 const int s13c = parameters.dbase.get<
int >(
"s13c");
220 const int s21c = parameters.dbase.get<
int >(
"s21c");
221 const int s22c = parameters.dbase.get<
int >(
"s22c");
222 const int s23c = parameters.dbase.get<
int >(
"s23c");
223 const int s31c = parameters.dbase.get<
int >(
"s31c");
224 const int s32c = parameters.dbase.get<
int >(
"s32c");
225 const int s33c = parameters.dbase.get<
int >(
"s33c");
227 const int pc = parameters.dbase.get<
int >(
"pc");
229 bool assignStress = s11c >=0 ;
235 real cs = sqrt(
mu/
rho );
237 std::vector<real> & twd = parameters.dbase.get<std::vector<real> >(
"travelingWaveData");
238 const int np = int(twd[0]);
239 const int ns = int(twd[1]);
243 printF(
"\n\n **************** FIX ME: travelingWave: finish me for HEMP **********\n\n");
251 if(
mg.numberOfDimensions()==2 )
256 real x0 =
X(i1,i2,i3,0);
257 real y0 =
X(i1,i2,i3,1);
261 real s11=0., s12=0., s22=0.;
268 for (
int n=0;
n<
np;
n++ )
270 real ap = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
273 real xi = k1*(x0-xa) + k2*(y0-ya) - cp*t;
274 real sinf = sin(freq*xi), cosf=cos(freq*xi);
278 v1 += ap*cp*k1*freq*cosf;
279 v2 += ap*cp*k2*freq*cosf;
280 s11+= -ap*(
lambda+2.*
mu*k1*k1 )*freq*cosf;
281 s12+= -ap*( 2.*
mu*k1*k2 )*freq*cosf;
282 s22+= -ap*(
lambda+2.*
mu*k2*k2 )*freq*cosf;
285 for (
int n=0;
n<
ns;
n++ )
287 real as = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
289 real xi = k1*(x0-xa) + k2*(y0-ya) - cs*t;
291 real sinf = sin(freq*xi), cosf=cos(freq*xi);
296 v1 += -as*cs*k2*freq*cosf;
297 v2 += +as*cs*k1*freq*cosf;
298 s11+= -as*(-2.*
mu*k1*k2 )*freq*cosf;
299 s12+= -as*(
mu*(k1*k1-k2*k2) )*freq*cosf;
300 s22+= -as*( 2.*
mu*k1*k2 )*freq*cosf;
315 state0Local(i1,i2,i3,1) = 1.0;
317 state0Local(i1,i2,i3,0) = state0Local(i1,i2,i3,1)*det(i1,i2,i3)*
mg.gridSpacing(axis1)*
mg.gridSpacing(axis2);
319 i1 > I1Base && i1 < I1Bound &&
320 i2 > I2Base && i2 < I2Bound )
322 real area = 0.25*(det(i1,i2,i3)+det(i1+1,i2,i3)+det(i1+1,i2+1,i3)+det(i1,i2+1,i3));
323 state0Local(i1,i2,i3,0) = area*
mg.gridSpacing(axis1)*
mg.gridSpacing(axis2);
333 if( assignVelocities )
335 U(i1,i2,i3,v1c) = v1;
336 U(i1,i2,i3,v2c) = v2;
340 U(i1,i2,i3,s11c) =s11;
341 U(i1,i2,i3,s12c) =s12;
342 U(i1,i2,i3,s21c) =s12;
343 U(i1,i2,i3,s22c) =s22;
346 real press = -(
lambda+2.0*
mu/3.0)*div;
347 U(i1,i2,i3,pc) = press;
348 U(i1,i2,i3,s11c) += press;
349 U(i1,i2,i3,s22c) += press;
358 ERR(i1,i2,i3,u1c) =
U(i1,i2,i3,u1c) - u1;
359 ERR(i1,i2,i3,u2c) =
U(i1,i2,i3,u2c) - u2;
360 ERR(i1,i2,i3,uc) =
U(i1,i2,i3,uc) - x0;
361 ERR(i1,i2,i3,vc) =
U(i1,i2,i3,vc) - y0;
365 ERR(i1,i2,i3,uc) =
U(i1,i2,i3,uc) - u1;
366 ERR(i1,i2,i3,vc) =
U(i1,i2,i3,vc) - u2;
369 if( assignVelocities )
371 ERR(i1,i2,i3,v1c) =
U(i1,i2,i3,v1c) - v1;
372 ERR(i1,i2,i3,v2c) =
U(i1,i2,i3,v2c) - v2;
376 ERR(i1,i2,i3,s11c) =
U(i1,i2,i3,s11c) -s11;
377 ERR(i1,i2,i3,s12c) =
U(i1,i2,i3,s12c) -s12;
378 ERR(i1,i2,i3,s21c) =
U(i1,i2,i3,s21c) -s12;
379 ERR(i1,i2,i3,s22c) =
U(i1,i2,i3,s22c) -s22;
382 real press = -(
lambda+2.0*
mu/3.0)*div;
383 ERR(i1,i2,i3,s11c) -= press;
384 ERR(i1,i2,i3,s22c) -= press;
385 ERR(i1,i2,i3,pc) =
U(i1,i2,i3,pc) - press;
411 #beginMacro getRayleighWaveSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
413 const int v1c = parameters.dbase.get<
int >(
"v1c");
414 const int v2c = parameters.dbase.get<
int >(
"v2c");
415 const int v3c = parameters.dbase.get<
int >(
"v3c");
417 bool assignVelocities= v1c>=0 ;
418 const int s11c = parameters.dbase.get<
int >(
"s11c");
419 const int s12c = parameters.dbase.get<
int >(
"s12c");
420 const int s13c = parameters.dbase.get<
int >(
"s13c");
421 const int s21c = parameters.dbase.get<
int >(
"s21c");
422 const int s22c = parameters.dbase.get<
int >(
"s22c");
423 const int s23c = parameters.dbase.get<
int >(
"s23c");
424 const int s31c = parameters.dbase.get<
int >(
"s31c");
425 const int s32c = parameters.dbase.get<
int >(
"s32c");
426 const int s33c = parameters.dbase.get<
int >(
"s33c");
428 const int pc = parameters.dbase.get<
int >(
"pc");
430 bool assignStress = s11c >=0 ;
433 real cs = sqrt(
mu/
rho );
435 std::vector<real> &
data = parameters.dbase.get<std::vector<real> >(
"RayleighWaveData");
436 const int nk = int(data[0]);
437 const real
cr = data[1];
445 printF(
"\n\n **************** FIX ME: RayelighWave: finish me for HEMP **********\n\n");
449 real
cb1 = sqrt(1.-SQR(cr/cp));
450 real
cb2 = sqrt(1.-SQR(cr/cs));
452 real
c1 = .5*SQR(cr/cs)-1.;
456 printF(
"**** RayleighWave: ySurf=%8.2e, cr=%8.2e, period=%8.2e, t=%8.2e *********\n",ySurf,cr,period,t);
458 for(
int n=0;
n<
nk;
n++ )
463 printF(
" k%i = %e, a%i=%e, b%i=%e\n",
n,k,
n,a,
n,b);
472 if(
mg.numberOfDimensions()==2 )
477 real x0 =
X(i1,i2,i3,0)-
xShift;
478 real y0 =
X(i1,i2,i3,1)-
ySurf;
482 real s11=0., s12=0., s22=0.;
486 for(
int n=0;
n<
nk;
n++ )
488 real
k = twoPi*data[m++]/
period;
491 real b=data[m++]*
scale;
496 real eb1 = exp(b1*(y0)), eb2 = exp(b2*(y0));
498 real ct = cos(k*(x0-cr*t));
499 real st = sin(k*(x0-cr*t));
501 u1 += ( eb1 + c1*eb2 )*( a*ct+b*st);
502 u2 += ( cb1*eb1 + (c1/
cb2)*eb2 )*( a*st-b*ct);
503 if( assignVelocities )
505 v1 += ( eb1 + c1*eb2 )*( k*cr*( a*st-b*ct) );
506 v2 +=-( cb1*eb1 + (c1/
cb2)*eb2 )*( k*cr*( a*ct+b*st) );
510 real u1x = ( eb1 + c1*eb2 )*( k*(-a*st+b*ct) );
511 real u1y = ( b1*eb1 + b2*c1*eb2 )*( ( a*ct+b*st) );
513 real u2x = ( cb1*eb1 + (c1/
cb2)*eb2 )*( k*(a*ct+b*st) );
514 real u2y = ( b1*cb1*eb1 + b2*(c1/
cb2)*eb2 )*( (a*st-b*ct) );
518 s12 +=
mu*( u1y+u2x );
532 if( assignVelocities )
534 U(i1,i2,i3,v1c) = v1;
535 U(i1,i2,i3,v2c) = v2;
539 U(i1,i2,i3,s11c) =s11;
540 U(i1,i2,i3,s12c) =s12;
541 U(i1,i2,i3,s21c) =s12;
542 U(i1,i2,i3,s22c) =s22;
549 ERR(i1,i2,i3,u1c) =
U(i1,i2,i3,u1c) - u1;
550 ERR(i1,i2,i3,u2c) =
U(i1,i2,i3,u2c) - u2;
551 ERR(i1,i2,i3,
uc) =
U(i1,i2,i3,
uc) - x0;
552 ERR(i1,i2,i3,
vc) =
U(i1,i2,i3,
vc) - y0;
556 ERR(i1,i2,i3,
uc) =
U(i1,i2,i3,
uc) - u1;
557 ERR(i1,i2,i3,
vc) =
U(i1,i2,i3,
vc) - u2;
560 if( assignVelocities )
562 ERR(i1,i2,i3,v1c) =
U(i1,i2,i3,v1c) - v1;
563 ERR(i1,i2,i3,v2c) =
U(i1,i2,i3,v2c) - v2;
567 ERR(i1,i2,i3,s11c) =
U(i1,i2,i3,s11c) -s11;
568 ERR(i1,i2,i3,s12c) =
U(i1,i2,i3,s12c) -s12;
569 ERR(i1,i2,i3,s21c) =
U(i1,i2,i3,s21c) -s12;
570 ERR(i1,i2,i3,s22c) =
U(i1,i2,i3,s22c) -s22;
580 OV_ABORT(
"RayleighWave:ERROR: finish me for 3D");
589 #defineMacro fg(x,z) ( .5*( cg1*(x)*( 1.+ (z)*( 7./p + (z)*( 21./(2.*p-1.) + (z)*( 35./(3.*p-2.) + (z)*( 35./(4.*p-3.) \
590 + (z)*( 21./(5.*p-4.) + (z)*( 7./(6.*p-5.) + z/(7.*p-6.) ))))))) ) )
597 #beginMacro getF( x, f, fPrime )
602 f = cp*fg(xx,z) - .5*(
a/p)*xp1*xx;
611 #beginMacro getG( x, g, gPrime )
616 real
xp1 = pow(xx,p-1);
618 g = cp*fg(xx,z) + .5*(
a/p)*xp1*xx ;
619 gPrime = +.5*( -cg1*pow( 1. + z , 7.) -
a*xp1/
cp );
625 real
xp1 = pow(xx,p-1);
627 g = -(
a/p)*xp1*xx - cp*fg(xx,z) + .5*(
a/p)*xp1*xx;
628 gPrime = -
a*xp1/cp + .5*( -cg1*pow( 1. + z , 7.) +
a*xp1/
cp );
638 #beginMacro getPistonMotionSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
640 const int v1c = parameters.dbase.get<
int >(
"v1c");
641 const int v2c = parameters.dbase.get<
int >(
"v2c");
642 const int v3c = parameters.dbase.get<
int >(
"v3c");
644 bool assignVelocities= v1c>=0 ;
645 const int s11c = parameters.dbase.get<
int >(
"s11c");
646 const int s12c = parameters.dbase.get<
int >(
"s12c");
647 const int s13c = parameters.dbase.get<
int >(
"s13c");
648 const int s21c = parameters.dbase.get<
int >(
"s21c");
649 const int s22c = parameters.dbase.get<
int >(
"s22c");
650 const int s23c = parameters.dbase.get<
int >(
"s23c");
651 const int s31c = parameters.dbase.get<
int >(
"s31c");
652 const int s32c = parameters.dbase.get<
int >(
"s32c");
653 const int s33c = parameters.dbase.get<
int >(
"s33c");
655 const int pc = parameters.dbase.get<
int >(
"pc");
657 bool assignStress = s11c >=0 ;
660 const real cs = sqrt(
mu/
rho );
662 std::vector<real> & data = parameters.dbase.get<std::vector<real> >(
"pistonMotionData");
665 const real
a =data[m++];
666 const real p =data[m++];
667 const real rhog =data[m++];
668 const real pg =data[m++];
669 const real gamma=data[m++];
670 real angle =data[m++];
672 const real
a0 = sqrt( gamma*pg/rhog);
676 printF(
"\n\n **************** FIX ME: getPistonMotionSolution: finish me for HEMP **********\n\n");
685 printP(
"**** getPistonMotion: a=%8.2e, p=%8.2e, Gas: rho=%8.2e, p=%8.2e gamma=%8.2e angle=%5.2f(degrees)*******\n",
686 a,p,rhog,pg,gamma,angle);
689 angle = angle*Pi/180.;
690 const real cosa = cos(angle), sina=sin(angle);
692 const real cg1 = pg/(
rho*cp*
cp);
693 const real cg2 = (-
a)*(gamma-1.)/(2.*
a0);
695 assert( fabs(gamma-1.4) < REAL_EPSILON*100. );
698 if(
mg.numberOfDimensions()==2 )
703 real x0 =
X(i1,i2,i3,0);
704 real y0 =
X(i1,i2,i3,1);
708 real s11=0., s12=0., s22=0.;
735 real xa = x0*cosa + y0*sina;
741 real fm, fmPrime, gp, gpPrime;
743 getF( xm, fm, fmPrime );
744 getG( xp, gp, gpPrime );
747 real va = cp*( - fmPrime + gpPrime );
748 real uap=fmPrime + gpPrime;
766 real u1x = uap*cosa*cosa;
767 real u1y = uap*cosa*sina;
769 real u2x = uap*sina*cosa;
770 real u2y = uap*sina*sina;
774 s12 =
mu*( u1y+u2x );
784 if( assignVelocities )
786 U(i1,i2,i3,v1c) = v1;
787 U(i1,i2,i3,v2c) = v2;
791 U(i1,i2,i3,s11c) =s11;
792 U(i1,i2,i3,s12c) =s12;
793 U(i1,i2,i3,s21c) =s12;
794 U(i1,i2,i3,s22c) =s22;
799 ERR(i1,i2,i3,uc) =
U(i1,i2,i3,uc) - u1;
800 ERR(i1,i2,i3,vc) =
U(i1,i2,i3,vc) - u2;
802 if( assignVelocities )
804 ERR(i1,i2,i3,v1c) =
U(i1,i2,i3,v1c) - v1;
805 ERR(i1,i2,i3,v2c) =
U(i1,i2,i3,v2c) - v2;
809 ERR(i1,i2,i3,s11c) =
U(i1,i2,i3,s11c) -s11;
810 ERR(i1,i2,i3,s12c) =
U(i1,i2,i3,s12c) -s12;
811 ERR(i1,i2,i3,s21c) =
U(i1,i2,i3,s21c) -s12;
812 ERR(i1,i2,i3,s22c) =
U(i1,i2,i3,s22c) -s22;
822 OV_ABORT(
"getPistonMotion:ERROR: finish me for 3D");
839 #define exTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(-ky/(eps*cc))
840 #define eyTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*( kx/(eps*cc))
841 #define hzTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))
843 #define exLaplacianTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(+ky*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
844 #define eyLaplacianTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(-kx*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
845 #define hzLaplacianTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*( -(twoPi*twoPi*(kx*kx+ky*ky) ) )
852 #define hzGaussianPulse(xi) exp(-betaGaussianPlaneWave*((xi)*(xi)))
853 #define exGaussianPulse(xi) hzGaussianPulse(xi)*(-ky/(eps*cc))
854 #define eyGaussianPulse(xi) hzGaussianPulse(xi)*( kx/(eps*cc))
856 #define hzLaplacianGaussianPulse(xi) ((4.*betaGaussianPlaneWave*betaGaussianPlaneWave*(kx*kx+ky*ky))*xi*xi-\
857 (2.*betaGaussianPlaneWave*(kx*kx+ky*ky)))*exp(-betaGaussianPlaneWave*((xi)*(xi)))
858 #define exLaplacianGaussianPulse(xi) hzLaplacianGaussianPulse(xi,t)*(-ky/(eps*cc))
859 #define eyLaplacianGaussianPulse(xi) hzLaplacianGaussianPulse(xi,t)*( kx/(eps*cc))
880 #define exTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)+kz*(z)-cc*(t)))*(-ky/(eps*cc))
881 #define eyTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)+kz*(z)-cc*(t)))*( kx/(eps*cc))
882 #define ezTrue3d(x,y,z,t) 0
884 #define hxTrue3d(x,y,z,t) 0
885 #define hyTrue3d(x,y,z,t) 0
886 #define hzTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)+kz*(z)-cc*(t)))
888 #define exLaplacianTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(+ky*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
889 #define eyLaplacianTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(-kx*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
890 #define ezLaplacianTrue3d(x,y,z,t) 0
892 #define hxLaplacianTrue3d(x,y,z,t) 0
893 #define hyLaplacianTrue3d(x,y,z,t) 0
894 #define hzLaplacianTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*( -(twoPi*twoPi*(kx*kx+ky*ky) ) )
898 #defineMacro PMIex(x,y,t) (aa1*cos(twoPi*(kx*(x)+ky*(y)-w*(t))) + r*bb1*cos(twoPi*(-kx*(x)+ky*(y)-w*(t))))
899 #defineMacro PMIey(x,y,t) (aa2*cos(twoPi*(kx*(x)+ky*(y)-w*(t))) + r*bb2*cos(twoPi*(-kx*(x)+ky*(y)-w*(t))))
900 #defineMacro PMIhz(x,y,t) ((-twoPi/w)*( (a1*ky-a2*kx)*sin(twoPi*(kx*(x)+ky*(y)-w*(t))) \
901 + (b1*ky+b2*kx)*r*sin(twoPi*(-kx*(x)+ky*(y)-w*(t)))) )
903 #defineMacro PMITex(x,y,t) (tau*dd1*cos(twoPi*(kappax*(x)+kappay*(y)-w*(t))))
904 #defineMacro PMITey(x,y,t) (tau*dd2*cos(twoPi*(kappax*(x)+kappay*(y)-w*(t))))
905 #defineMacro PMIThz(x,y,t) ((-twoPi*tau/w)*(d1*kappay-d2*kappax)*sin(twoPi*(kappax*(x)+kappay*(y)-w*(t))))
908 #beginMacro setPlaneMaterialInterfaceMacro(OPTION,J1,J2,J3)
909 if( numberOfDimensions==2 )
914 const real c1=cGrid(0);
915 const real
c2=cGrid(1);
917 const real kNorm = sqrt(kx*kx+ky*ky);
918 const real a1=-ky/kNorm,
a2= kx/kNorm;
919 const real b1=-ky/kNorm, b2=-kx/kNorm;
920 const real w = c1*kNorm;
922 real kappax,kappay,kappaNorm;
927 kappax=sqrt( (kx*kx+ky*ky)/(cr*cr) - kappay*kappay );
929 kappaNorm=sqrt(kappax*kappax+kappay*kappay);
930 const real d1=-kappay/kappaNorm;
931 const real d2= kappax/kappaNorm;
933 const real cosTheta1=kx/kNorm;
934 const real cosTheta2=kappax/kappaNorm;
937 const real
r = (c1*cosTheta1-c2*cosTheta2)/(c1*cosTheta1+c2*cosTheta2);
938 const real tau = (2.*c2*cosTheta1)/(c1*cosTheta1+c2*cosTheta2);
941 const real x0=x0PlaneMaterialInterface[0], y0=x0PlaneMaterialInterface[1];
943 real an0=normalPlaneMaterialInterface[0], an1=normalPlaneMaterialInterface[1];
944 real aNorm = sqrt(an0*an0+an1*an1);
945 an0/= aNorm; an1/=aNorm;
946 real aa1=an0*a1-an1*
a2, aa2= an1*a1+an0*
a2;
947 real bb1=an0*b1-an1*b2, bb2= an1*b1+an0*b2;
948 real dd1=an0*d1-an1*d2, dd2= an1*d1+an0*d2;
954 real xx1 =
XEP(i1,i2,i3,0) -x0;
955 real yy1 =
XEP(i1,i2,i3,1) -y0;
960 real u1 = PMIex(x,y,t);
961 real u2 = PMIey(x,y,t);
962 real u3 = PMIhz(x,y,t);
964 #If #OPTION eq "initialCondition"
969 UMEX(i1,i2,i3)= PMIex(x,y,tm);
970 UMEY(i1,i2,i3)= PMIey(x,y,tm);
971 UMHZ(i1,i2,i3)= PMIhz(x,y,tm);
972 #Elif #OPTION eq "error"
976 #Elif #OPTION eq "boundaryCondition"
988 real xx1 =
XEP(i1,i2,i3,0) -x0;
989 real yy1 =
XEP(i1,i2,i3,1) -y0;
994 real u1 = PMITex(x,y,t);
995 real u2 = PMITey(x,y,t);
996 real u3 = PMIThz(x,y,t);
998 #If #OPTION eq "initialCondition"
1003 UMEX(i1,i2,i3)= PMITex(x,y,tm);
1004 UMEY(i1,i2,i3)= PMITey(x,y,tm);
1005 UMHZ(i1,i2,i3)= PMIThz(x,y,tm);
1006 #Elif #OPTION eq "error"
1010 #Elif #OPTION eq "boundaryCondition"
1027 #beginMacro getGaussianIntegralSolution(OPTION,VEX,VEY,VHZ,t)
1033 const int nsources=1;
1034 double xs[nsources], ys[nsources], tau[nsources], var[nsources], amp[nsources];
1048 double x=
X(i1,i2,i3,0);
1049 double y=
X(i1,i2,i3,1);
1051 exmax(wt,wx,wy,nsources,xs[0],ys[0],tau[0],var[0],amp[0],period,x,y,time);
1053 #If #OPTION eq "solution"
1057 #Elif #OPTION eq "error"
1058 ERREX(i1,i2,i3) = VEX(i1,i2,i3) - wy;
1059 ERREY(i1,i2,i3) = VEY(i1,i2,i3) + wx;
1060 ERRHZ(i1,i2,i3) = VHZ(i1,i2,i3) - wt;
1063 U(i1,i2,i3,
ex) = wy;
1064 U(i1,i2,i3,
ey) =-wx;
1065 U(i1,i2,i3,hz) = wt;