00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "Field/BCond.h"
00028 #include "Field/BareField.h"
00029 #include "Index/NDIndex.h"
00030 #include "Index/Index.h"
00031 #include "Field/GuardCellSizes.h"
00032 #include "Field/BrickIterator.h"
00033 #include "Field/BrickExpression.h"
00034 #include "Meshes/Centering.h"
00035 #include "Meshes/CartesianCentering.h"
00036 #include "Utility/IpplInfo.h"
00037 #include "Utility/PAssert.h"
00038 #include "AppTypes/AppTypeTraits.h"
00039 #include "Profile/Profiler.h"
00040
00041 #ifdef IPPL_USE_STANDARD_HEADERS
00042 #include <iostream>
00043 #include <typeinfo>
00044 #include <vector>
00045 using namespace std;
00046 #else
00047 #include <vector.h>
00048 #include <iostream.h>
00049 #include <typeinfo.h>
00050 #endif
00051
00053
00054 template<class T, unsigned D, class M, class C>
00055 int BCondBase<T,D,M,C>::allComponents = -9999;
00056
00058
00059
00060
00061
00062 #define COMPONENT_APPLY_BUILTIN(OP,T) \
00063 inline void PETE_apply(const OP<T>&, T&, const T&) \
00064 { \
00065 ERRORMSG("Component boundary condition on a scalar (T)." << endl); \
00066 Ippl::abort(); \
00067 }
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077 template<class T, unsigned int D, class M, class C>
00078 BCondBase<T,D,M,C>::BCondBase(unsigned int face, int i, int j)
00079 : m_face(face), m_changePhysical(false)
00080 {
00081
00082
00083
00084 if (j != BCondBase<T,D,M,C>::allComponents) {
00085 if (i == BCondBase<T,D,M,C>::allComponents)
00086 ERRORMSG("BCondBase(): component 2 specified, component 1 not."
00087 << endl);
00088
00089
00090
00091
00092
00093 if (getTensorOrder(get_tag(T())) == IPPL_TENSOR) {
00094 m_component = i + j*D;
00095 } else if (getTensorOrder(get_tag(T())) == IPPL_SYMTENSOR) {
00096 int lo = i < j ? i : j;
00097 int hi = i > j ? i : j;
00098 m_component = ((hi+1)*hi/2) + lo;
00099 } else if (getTensorOrder(get_tag(T())) == IPPL_ANTISYMTENSOR) {
00100 PAssert(i > j);
00101 m_component = ((i-1)*i/2) + j;
00102 } else {
00103 ERRORMSG(
00104 "BCondBase(): something other than [Sym,AntiSym]Tenzor specified"
00105 << " two component indices; not implemented." << endl);
00106 }
00107
00108 } else {
00109
00110
00111
00112
00113 m_component = i;
00114 }
00115 }
00116
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 template<class T, unsigned int D, class M, class C>
00129 void BCondBase<T,D,M,C>::write(ostream& out) const
00130 {
00131 out << "BCondBase" << ", Face=" << m_face;
00132 }
00133
00134 template<class T, unsigned int D, class M, class C>
00135 void PeriodicFace<T,D,M,C>::write(ostream& out) const
00136 {
00137 out << "PeriodicFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00138 }
00139
00140 template<class T, unsigned int D, class M, class C>
00141 void ParallelPeriodicFace<T,D,M,C>::write(ostream& out) const
00142 {
00143 out << "ParallelPeriodicFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00144 }
00145
00146 template<class T, unsigned int D, class M, class C>
00147 void NegReflectFace<T,D,M,C>::write(ostream& out) const
00148 {
00149 out << "NegReflectFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00150 }
00151
00152 template<class T, unsigned int D, class M, class C>
00153 void NegReflectAndZeroFace<T,D,M,C>::write(ostream& out) const
00154 {
00155 out << "NegReflectAndZeroFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00156 }
00157
00158 template<class T, unsigned int D, class M, class C>
00159 void PosReflectFace<T,D,M,C>::write(ostream& out) const
00160 {
00161 out << "PosReflectFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00162 }
00163
00164 template<class T, unsigned int D, class M, class C>
00165 void ZeroFace<T,D,M,C>::write(ostream& out) const
00166 {
00167 out << "ZeroFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00168 }
00169
00170 template<class T, unsigned int D, class M, class C>
00171 void ZeroGuardsAndZeroFace<T,D,M,C>::write(ostream& out) const
00172 {
00173 out << "ZeroGuardsAndZeroFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00174 }
00175
00176 template<class T, unsigned int D, class M, class C>
00177 void ConstantFace<T,D,M,C>::write(ostream& out) const
00178 {
00179 out << "ConstantFace"
00180 << ", Face=" << BCondBase<T,D,M,C>::m_face
00181 << ", Constant=" << this->Offset
00182 << endl;
00183 }
00184
00185 template<class T, unsigned int D, class M, class C>
00186 void EurekaFace<T,D,M,C>::write(ostream& out) const
00187 {
00188 out << "EurekaFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00189 }
00190
00191 template<class T, unsigned int D, class M, class C>
00192 void FunctionFace<T,D,M,C>::write(ostream& out) const
00193 {
00194 out << "FunctionFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00195 }
00196
00197 template<class T, unsigned int D, class M, class C>
00198 void ComponentFunctionFace<T,D,M,C>::write(ostream& out) const
00199 {
00200 out << "ComponentFunctionFace" << ", Face=" << BCondBase<T,D,M,C>::m_face;
00201 }
00202
00203 template<class T, unsigned D, class M, class C>
00204 void
00205 ExtrapolateFace<T,D,M,C>::write(ostream& o) const
00206 {
00207 TAU_TYPE_STRING(taustr, CT(*this) + " void (ostream )" );
00208 TAU_PROFILE("ExtrapolateFace::write()", taustr, TAU_FIELD | TAU_IO);
00209 o << "ExtrapolateFace, Face=" << BCondBase<T,D,M,C>::m_face
00210 << ", Offset=" << Offset << ", Slope=" << Slope;
00211 }
00212
00213 template<class T, unsigned D, class M, class C>
00214 void
00215 ExtrapolateAndZeroFace<T,D,M,C>::write(ostream& o) const
00216 {
00217 TAU_TYPE_STRING(taustr, CT(*this) + " void (ostream )" );
00218 TAU_PROFILE("ExtrapolateAndZeroFace::write()", taustr, TAU_FIELD | TAU_IO);
00219 o << "ExtrapolateAndZeroFace, Face=" << BCondBase<T,D,M,C>::m_face
00220 << ", Offset=" << Offset << ", Slope=" << Slope;
00221 }
00222
00223 template<class T, unsigned D, class M, class C>
00224 void
00225 LinearExtrapolateFace<T,D,M,C>::write(ostream& o) const
00226 {
00227 TAU_TYPE_STRING(taustr, CT(*this) + " void (ostream )" );
00228 TAU_PROFILE("LinearExtrapolateFace::write()", taustr, TAU_FIELD | TAU_IO);
00229 o << "LinearExtrapolateFace, Face=" << BCondBase<T,D,M,C>::m_face;
00230 }
00231
00232 template<class T, unsigned D, class M, class C>
00233 void
00234 ComponentLinearExtrapolateFace<T,D,M,C>::write(ostream& o) const
00235 {
00236 TAU_TYPE_STRING(taustr, CT(*this) + " void (ostream )" );
00237 TAU_PROFILE("ComponentLinearExtrapolateFace::write()",
00238 taustr, TAU_FIELD | TAU_IO);
00239 o << "ComponentLinearExtrapolateFace, Face=" << BCondBase<T,D,M,C>::m_face;
00240 }
00241
00243
00244 template<class T, unsigned D, class M, class C>
00245 void
00246 BConds<T,D,M,C>::write(ostream& o) const
00247 {
00248 TAU_TYPE_STRING(taustr, CT(*this) + " void (ostream )" );
00249 TAU_PROFILE("BConds::write()", taustr, TAU_FIELD | TAU_IO);
00250
00251 o << "BConds:(" << endl;
00252 const_iterator p=this->begin();
00253 while (p!=this->end())
00254 {
00255 (*p).second->write(o);
00256 ++p;
00257 if (p!=this->end())
00258 o << " , " << endl;
00259 else
00260 o << endl << ")" << endl << endl;
00261 }
00262 }
00263
00265
00266 template<class T, unsigned D, class M, class C>
00267 void
00268 BConds<T,D,M,C>::apply( Field<T,D,M,C>& a )
00269 {
00270 TAU_TYPE_STRING(taustr, "void (" + CT(a) + " )" );
00271 TAU_PROFILE("BConds::apply()", taustr, TAU_FIELD | TAU_ASSIGN);
00272 for (iterator p=this->begin(); p!=this->end(); ++p)
00273 (*p).second->apply(a);
00274 }
00275
00276 template<class T, unsigned D, class M, class C>
00277 bool
00278 BConds<T,D,M,C>::changesPhysicalCells() const
00279 {
00280 for (const_iterator p=this->begin(); p!=this->end(); ++p)
00281 if ((*p).second->changesPhysicalCells())
00282 return true;
00283 return false;
00284 }
00285
00286
00287
00288
00289
00290
00291 template<class T, unsigned D, class M, class C>
00292 PeriodicFace<T,D,M,C>::PeriodicFace(unsigned f, int i, int j)
00293 : BCondBase<T,D,M,C>(f,i,j)
00294 {
00295 TAU_TYPE_STRING(taustr, CT(*this) + " void (unsigned, int, int)" );
00296 TAU_PROFILE("PeriodicFace::PeriodicFace()", taustr, TAU_FIELD);
00297 }
00298
00299 template<class T, unsigned D, class M, class C>
00300 ExtrapolateFace<T,D,M,C>::ExtrapolateFace(unsigned f, T o, T s, int i, int j)
00301 : BCondBase<T,D,M,C>(f,i,j), Offset(o), Slope(s)
00302 {
00303 TAU_TYPE_STRING(taustr, "void (unsigned, " + CT(o) + ", " + CT(s) + ", int, int )" );
00304 TAU_PROFILE("ExtrapolateFace::ExtrapolateFace()", taustr, TAU_FIELD);
00305 }
00306
00307 template<class T, unsigned D, class M, class C>
00308 ExtrapolateAndZeroFace<T,D,M,C>::
00309 ExtrapolateAndZeroFace(unsigned f, T o, T s, int i, int j)
00310 : BCondBase<T,D,M,C>(f,i,j), Offset(o), Slope(s)
00311 {
00312 TAU_TYPE_STRING(taustr,
00313 "void (unsigned, " + CT(o) + ", " + CT(s) + ", int, int )" );
00314 TAU_PROFILE("ExtrapolateAndZeroFace::ExtrapolateAndZeroFace()",
00315 taustr, TAU_FIELD);
00316 BCondBase<T,D,M,C>::m_changePhysical = true;
00317 }
00318
00319 template<class T, unsigned D, class M, class C>
00320 FunctionFace<T,D,M,C>::
00321 FunctionFace(T (*func)(const T&), unsigned face)
00322 : BCondBase<T,D,M,C>(face), Func(func)
00323 {
00324 TAU_TYPE_STRING(taustr, "void (" + CT(func) + ", unsigned )" );
00325 TAU_PROFILE("FunctionFace::FunctionFace()", taustr, TAU_FIELD);
00326 }
00327
00328 template<class T, unsigned D, class M, class C>
00329 ComponentFunctionFace<T,D,M,C>::
00330 ComponentFunctionFace(typename ApplyToComponentType<T>::type
00331 (*func)( typename ApplyToComponentType<T>::type),
00332 unsigned face, int i, int j)
00333 : BCondBase<T,D,M,C>(face,i,j), Func(func)
00334 {
00335 TAU_TYPE_STRING(taustr, "void (" + CT(func) + ", unsigned, int, int )" );
00336 TAU_PROFILE("ComponentFunctionFace::ComponentFunctionFace()", taustr,
00337 TAU_FIELD | TAU_ASSIGN);
00338
00339
00340 if ( (j == BCondBase<T,D,M,C>::allComponents) &&
00341 (i == BCondBase<T,D,M,C>::allComponents) )
00342 ERRORMSG("ComponentFunctionFace(): allComponents specified; not allowed; "
00343 << "use FunctionFace class instead." << endl);
00344 }
00345
00346
00348
00349
00350
00351
00352
00353
00354 template<class T>
00355 struct OpPeriodic
00356 {
00357 #ifdef IPPL_PURIFY
00358 OpPeriodic() {}
00359 OpPeriodic(const OpPeriodic<T> &) {}
00360 OpPeriodic<T>& operator=(const OpPeriodic<T> &) { return *this; }
00361 #endif
00362 };
00363 template<class T>
00364 inline void PETE_apply(const OpPeriodic<T>& e, T& a, const T& b) { a = b; }
00365
00366
00367 template<class T>
00368 struct OpPeriodicComponent
00369 {
00370 OpPeriodicComponent(int c) : Component(c) {}
00371 int Component;
00372 };
00373
00374 template<class T>
00375 inline void PETE_apply(const OpPeriodicComponent<T>& e, T& a, const T& b)
00376 { a[e.Component] = b[e.Component]; }
00377
00378
00379
00380
00381
00382
00383
00384
00385 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,char)
00386 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,bool)
00387 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,int)
00388 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,unsigned)
00389 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,short)
00390 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,long)
00391 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,float)
00392 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,double)
00393 COMPONENT_APPLY_BUILTIN(OpPeriodicComponent,dcomplex)
00394
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407 template<class T, unsigned D, class M>
00408 void PeriodicFaceBCApply(PeriodicFace<T,D,M,Cell>& pf,
00409 Field<T,D,M,Cell>& A );
00410 template<class T, unsigned D, class M>
00411 void PeriodicFaceBCApply(PeriodicFace<T,D,M,Vert>& pf,
00412 Field<T,D,M,Vert>& A );
00413 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
00414 void PeriodicFaceBCApply(PeriodicFace<T,D,M,
00415 CartesianCentering<CE,D,NC> >& pf,
00416 Field<T,D,M,CartesianCentering<CE,D,NC> >& A );
00417
00418 template<class T, unsigned D, class M, class C>
00419 void PeriodicFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
00420 {
00421 TAU_TYPE_STRING(taustr, "void (" + CT(A) + " )" );
00422 TAU_PROFILE("PeriodicFace()", taustr, TAU_FIELD);
00423 PeriodicFaceBCApply(*this, A);
00424 }
00425
00426
00427
00428
00429
00430
00431 template<class T, unsigned D, class M>
00432 void PeriodicFaceBCApply(PeriodicFace<T,D,M,Cell>& pf,
00433 Field<T,D,M,Cell>& A )
00434 {
00435 TAU_TYPE_STRING(taustr, "void (" + CT(pf) + ", " + CT(A) + " )" );
00436 TAU_PROFILE("PeriodicFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
00437
00438
00439
00440
00441 const NDIndex<D>& domain( A.getDomain() );
00442 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
00443 unsigned d = pf.getFace()/2;
00444 int offset;
00445 if ( pf.getFace() & 1 )
00446 {
00447 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.leftGuard(d) );
00448 offset = -domain[d].length();
00449 }
00450 else
00451 {
00452 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
00453 offset = domain[d].length();
00454 }
00455
00456
00457 typename Field<T,D,M,Cell>::iterator_if fill_i;
00458 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
00459 {
00460
00461 LField<T,D> &fill = *(*fill_i).second;
00462 const NDIndex<D> &fill_alloc = fill.getAllocated();
00463 if ( slab.touches( fill_alloc ) )
00464 {
00465
00466 NDIndex<D> dest = slab.intersect( fill_alloc );
00467
00468
00469 NDIndex<D> src = dest;
00470 src[d] = src[d] + offset;
00471
00472
00473 typename Field<T,D,M,Cell>::iterator_if from_i;
00474 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
00475 {
00476
00477 LField<T,D> &from = *(*from_i).second;
00478 const NDIndex<D> &from_owned = from.getOwned();
00479 const NDIndex<D> &from_alloc = from.getAllocated();
00480
00481 if ( src.touches( from_owned ) )
00482 {
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495 fill.Uncompress();
00496
00497 NDIndex<D> from_it = src.intersect( from_alloc );
00498 NDIndex<D> fill_it = dest.plugBase( from_it );
00499
00500 typedef typename LField<T,D>::iterator LFI;
00501 LFI lhs = fill.begin(fill_it);
00502 LFI rhs = from.begin(from_it);
00503
00504
00505 if (pf.getComponent() == BCondBase<T,D,M,Cell>::allComponents) {
00506 BrickExpression<D,LFI,LFI,OpPeriodic<T> >
00507 (lhs,rhs,OpPeriodic<T>()).apply();
00508 } else {
00509 BrickExpression<D,LFI,LFI,OpPeriodicComponent<T> >
00510 (lhs,rhs,OpPeriodicComponent<T>(pf.getComponent())).apply();
00511 }
00512 }
00513 }
00514 }
00515 }
00516 }
00517
00518
00519
00520
00521
00522 template<class T, unsigned D, class M>
00523 void PeriodicFaceBCApply(PeriodicFace<T,D,M,Vert>& pf,
00524 Field<T,D,M,Vert>& A )
00525 {
00526 TAU_TYPE_STRING(taustr, "void (" + CT(pf) + ", " + CT(A) + " )" );
00527 TAU_PROFILE("PeriodicFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
00528
00529
00530
00531
00532
00533 const NDIndex<D>& domain( A.getDomain() );
00534 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
00535 unsigned d = pf.getFace()/2;
00536 int offset;
00537 if ( pf.getFace() & 1 )
00538 {
00539
00540 slab[d] = Index(domain[d].max(), domain[d].max() + A.rightGuard(d));
00541
00542
00543 offset = -domain[d].length() + 1;
00544 }
00545 else
00546 {
00547 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
00548
00549
00550 offset = domain[d].length() - 1;
00551 }
00552
00553
00554 typename Field<T,D,M,Vert>::iterator_if fill_i;
00555 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
00556 {
00557
00558 LField<T,D> &fill = *(*fill_i).second;
00559 const NDIndex<D> &fill_alloc = fill.getAllocated();
00560 if ( slab.touches( fill_alloc ) )
00561 {
00562
00563 NDIndex<D> dest = slab.intersect( fill_alloc );
00564
00565
00566 NDIndex<D> src = dest;
00567 src[d] = src[d] + offset;
00568
00569
00570 typename Field<T,D,M,Vert>::iterator_if from_i;
00571 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
00572 {
00573
00574 LField<T,D> &from = *(*from_i).second;
00575 const NDIndex<D> &from_owned = from.getOwned();
00576 const NDIndex<D> &from_alloc = from.getAllocated();
00577
00578 if ( src.touches( from_owned ) )
00579 {
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592 fill.Uncompress();
00593
00594 NDIndex<D> from_it = src.intersect( from_alloc );
00595 NDIndex<D> fill_it = dest.plugBase( from_it );
00596
00597 typedef typename LField<T,D>::iterator LFI;
00598 LFI lhs = fill.begin(fill_it);
00599 LFI rhs = from.begin(from_it);
00600
00601
00602 if (pf.getComponent() == BCondBase<T,D,M,Vert>::allComponents) {
00603 BrickExpression<D,LFI,LFI,OpPeriodic<T> >
00604 (lhs,rhs,OpPeriodic<T>()).apply();
00605 } else {
00606 BrickExpression<D,LFI,LFI,OpPeriodicComponent<T> >
00607 (lhs,rhs,OpPeriodicComponent<T>(pf.getComponent())).apply();
00608 }
00609 }
00610 }
00611 }
00612 }
00613 }
00614
00615
00616
00617
00618
00619
00620 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
00621 void PeriodicFaceBCApply(PeriodicFace<T,D,M,
00622 CartesianCentering<CE,D,NC> >& pf,
00623 Field<T,D,M,CartesianCentering<CE,D,NC> >& A )
00624 {
00625 TAU_TYPE_STRING(taustr, "void (" + CT(pf) + ", " + CT(A) + " )" );
00626 TAU_PROFILE("PeriodicFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
00627
00628
00629
00630
00631
00632 const NDIndex<D>& domain( A.getDomain() );
00633 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
00634 unsigned d = pf.getFace()/2;
00635 int offset;
00636 if ( pf.getFace() & 1 )
00637 {
00638
00639
00640 if (pf.getComponent() == BCondBase<T,D,M,CartesianCentering<CE,D,NC> >::
00641 allComponents) {
00642
00643 CenteringEnum centering0 = CE[0 + d*NC];
00644 for (int c=1; c<NC; c++) {
00645 if (CE[c + d*NC] != centering0)
00646 ERRORMSG("PeriodicFaceBCApply: BCond thinks all components have"
00647 << " same centering along direction " << d
00648 << ", but it isn't so." << endl);
00649 }
00650
00651 if (centering0 == CELL) {
00652 offset = -domain[d].length();
00653 } else {
00654
00655 slab[d] =
00656 Index( domain[d].max(), domain[d].max() + A.rightGuard(d));
00657 offset = -domain[d].length()+1;
00658 }
00659 } else {
00660
00661 if (CE[pf.getComponent() + d*NC] == CELL) {
00662 offset = -domain[d].length();
00663 } else {
00664 slab[d] =
00665 Index( domain[d].max(), domain[d].max() + A.rightGuard(d));
00666 offset = -domain[d].length()+1;
00667 }
00668 }
00669 }
00670 else
00671 {
00672 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
00673
00674
00675 if (pf.getComponent() == BCondBase<T,D,M,CartesianCentering<CE,D,NC> >::
00676 allComponents) {
00677
00678 CenteringEnum centering0 = CE[0 + d*NC];
00679 for (int c=1; c<NC; c++) {
00680 if (CE[c + d*NC] != centering0)
00681 ERRORMSG("PeriodicFaceBCApply: BCond thinks all components have"
00682 << " same centering along direction " << d
00683 << ", but it isn't so." << endl);
00684 }
00685
00686 if (centering0 == CELL) {
00687 offset = -domain[d].length();
00688 } else {
00689 offset = -domain[d].length() + 1;
00690 }
00691 } else {
00692
00693 if (CE[pf.getComponent() + d*NC] == CELL) {
00694 offset = domain[d].length();
00695 } else {
00696 offset = domain[d].length() - 1;
00697 }
00698 }
00699 }
00700
00701
00702 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if fill_i;
00703 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
00704 {
00705
00706 LField<T,D> &fill = *(*fill_i).second;
00707 const NDIndex<D> &fill_alloc = fill.getAllocated();
00708 if ( slab.touches( fill_alloc ) )
00709 {
00710
00711 NDIndex<D> dest = slab.intersect( fill_alloc );
00712
00713
00714 NDIndex<D> src = dest;
00715 src[d] = src[d] + offset;
00716
00717
00718 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if from_i;
00719 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
00720 {
00721
00722 LField<T,D> &from = *(*from_i).second;
00723 const NDIndex<D> &from_owned = from.getOwned();
00724 const NDIndex<D> &from_alloc = from.getAllocated();
00725
00726 if ( src.touches( from_owned ) )
00727 {
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740 fill.Uncompress();
00741
00742 NDIndex<D> from_it = src.intersect( from_alloc );
00743 NDIndex<D> fill_it = dest.plugBase( from_it );
00744
00745 typedef typename LField<T,D>::iterator LFI;
00746 LFI lhs = fill.begin(fill_it);
00747 LFI rhs = from.begin(from_it);
00748
00749
00750 if (pf.getComponent() == BCondBase<T,D,M,
00751 CartesianCentering<CE,D,NC> >::allComponents) {
00752 BrickExpression<D,LFI,LFI,OpPeriodic<T> >
00753 (lhs,rhs,OpPeriodic<T>()).apply();
00754 } else {
00755 BrickExpression<D,LFI,LFI,OpPeriodicComponent<T> >
00756 (lhs,rhs,OpPeriodicComponent<T>(pf.getComponent())).apply();
00757 }
00758 }
00759 }
00760 }
00761 }
00762 }
00763
00764
00765
00766
00767
00768
00769
00770 #ifdef PRINT_DEBUG
00771
00772 # include <iterator.h>
00773 #endif
00774
00775 template <class T, unsigned D, class M>
00776 inline void
00777 CalcParallelPeriodicDomain(const Field<T,D,M,Cell> &A,
00778 const ParallelPeriodicFace<T,D,M,Cell>& pf,
00779 NDIndex<D> &dest_slab,
00780 int &offset)
00781 {
00782
00783
00784
00785 unsigned d = pf.getFace()/2;
00786
00787 const NDIndex<D>& domain(A.getDomain());
00788
00789 if (pf.getFace() & 1)
00790 {
00791
00792
00793
00794
00795
00796 dest_slab[d] =
00797 Index(domain[d].max() + 1, domain[d].max() + A.leftGuard(d));
00798
00799
00800
00801
00802
00803 offset = -domain[d].length();
00804 }
00805 else
00806 {
00807
00808
00809
00810
00811 dest_slab[d] =
00812 Index(domain[d].min() - A.leftGuard(d), domain[d].min()-1);
00813
00814
00815
00816 offset = domain[d].length();
00817 }
00818 }
00819
00820
00821
00822
00823 template <class T, unsigned D, class M>
00824 inline void
00825 CalcParallelPeriodicDomain(const Field<T,D,M,Vert> &A,
00826 const ParallelPeriodicFace<T,D,M,Vert>& pf,
00827 NDIndex<D> &dest_slab,
00828 int &offset)
00829 {
00830
00831
00832
00833 const NDIndex<D>& domain(A.getDomain());
00834
00835 unsigned d = pf.getFace()/2;
00836
00837 if (pf.getFace() & 1)
00838 {
00839
00840
00841
00842
00843
00844
00845 dest_slab[d] =
00846 Index(domain[d].max(), domain[d].max() + A.rightGuard(d));
00847
00848
00849
00850
00851
00852 offset = -domain[d].length() + 1;
00853 }
00854 else
00855 {
00856
00857
00858
00859
00860 dest_slab[d] =
00861 Index(domain[d].min() - A.leftGuard(d), domain[d].min()-1);
00862
00863
00864
00865 offset = domain[d].length() - 1;
00866 }
00867 }
00868
00869
00870
00871 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
00872 inline void
00873 CalcParallelPeriodicDomain(const Field<T,D,M,CartesianCentering<CE,D,NC> >& A,
00874 const ParallelPeriodicFace<T,D,M,
00875 CartesianCentering<CE,D,NC> >& pf,
00876 NDIndex<D> &dest_slab,
00877 int &offset)
00878 {
00879 typedef BCondBase<T,D,M,CartesianCentering<CE,D,NC> > BCBase_t;
00880
00881
00882
00883
00884 const NDIndex<D>& domain(A.getDomain());
00885
00886 unsigned d = pf.getFace()/2;
00887
00888 if (pf.getFace() & 1)
00889 {
00890
00891
00892
00893
00894
00895
00896 if (pf.getComponent() == BCBase_t::allComponents)
00897 {
00898
00899
00900
00901 CenteringEnum centering0 = CE[0 + d*NC];
00902
00903 for (int c = 1; c < NC; c++)
00904 {
00905
00906 if (CE[c + d*NC] != centering0)
00907 ERRORMSG("ParallelPeriodicFaceBCApply:"
00908 << "BCond thinks all components have"
00909 << " same centering along direction " << d
00910 << ", but it isn't so." << endl);
00911 }
00912
00913
00914
00915
00916 if (centering0 == CELL) {
00917 offset = -domain[d].length();
00918 } else {
00919 dest_slab[d] =
00920 Index(domain[d].max(), domain[d].max() + A.leftGuard(d));
00921 offset = -domain[d].length() + 1;
00922 }
00923
00924 }
00925 else
00926 {
00927
00928
00929
00930 if (CE[pf.getComponent() + d*NC] == CELL)
00931 {
00932 offset = -domain[d].length();
00933 }
00934 else
00935 {
00936 dest_slab[d] =
00937 Index(domain[d].max(), domain[d].max() + A.leftGuard(d));
00938 offset = -domain[d].length() + 1;
00939 }
00940 }
00941 }
00942 else
00943 {
00944
00945
00946
00947
00948 dest_slab[d] =
00949 Index(domain[d].min() - A.leftGuard(d), domain[d].min()-1);
00950
00951
00952
00953 if (pf.getComponent() == BCBase_t::allComponents)
00954 {
00955
00956
00957
00958 CenteringEnum centering0 = CE[0 + d*NC];
00959
00960 for (int c = 1; c < NC; c++)
00961 {
00962 if (CE[c + d*NC] != centering0)
00963 ERRORMSG("ParallelPeriodicFaceBCApply:"
00964 << "BCond thinks all components have"
00965 << " same centering along direction " << d
00966 << ", but it isn't so." << endl);
00967 }
00968
00969
00970
00971
00972 if (centering0 == CELL) {
00973 offset = -domain[d].length();
00974 } else {
00975 offset = -domain[d].length() + 1;
00976 }
00977
00978 }
00979 else
00980 {
00981
00982
00983
00984 if (CE[pf.getComponent() + d*NC] == CELL)
00985 {
00986 offset = domain[d].length();
00987 }
00988 else
00989 {
00990 offset = domain[d].length() - 1;
00991 }
00992 }
00993 }
00994 }
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006 template<class T, unsigned D, class M, class C>
01007 void ParallelPeriodicFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
01008 {
01009 TAU_TYPE_STRING(taustr, "void (" + CT(A) + " )" );
01010 TAU_PROFILE("ParallelPeriodicFace::apply()", taustr, TAU_FIELD);
01011
01012 TAU_PROFILE_TIMER(sendtimer,
01013 " ParallelPeriodicFaceBCApply-send",
01014 taustr, TAU_FIELD);
01015 TAU_PROFILE_TIMER(sendcommtimer,
01016 " ParallelPeriodicFaceBCApply-send-comm",
01017 taustr, TAU_FIELD);
01018 TAU_PROFILE_TIMER(findrectimer,
01019 " ParallelPeriodicFaceBCApply-findreceive",
01020 taustr, TAU_FIELD);
01021 TAU_PROFILE_TIMER(localstimer,
01022 " ParallelPeriodicFaceBCApply-locals",
01023 taustr, TAU_FIELD);
01024 TAU_PROFILE_TIMER(rectimer,
01025 " ParallelPeriodicFaceBCApply-receive",
01026 taustr, TAU_FIELD);
01027 TAU_PROFILE_TIMER(reccommtimer,
01028 " ParallelPeriodicFaceBCApply-receive-comm",
01029 taustr, TAU_FIELD);
01030 TAU_PROFILE_TIMER(localsexprtimer,
01031 " ParallelPeriodicFaceBCApply-locals-expression",
01032 taustr, TAU_FIELD);
01033
01034 #ifdef PRINT_DEBUG
01035 Inform msg("PPeriodicBC", INFORM_ALL_NODES);
01036 #endif
01037
01038
01039
01040
01041 typedef T Element_t;
01042 typedef NDIndex<D> Domain_t;
01043 typedef LField<T,D> LField_t;
01044 typedef typename LField_t::iterator LFI_t;
01045 typedef Field<T,D,M,C> Field_t;
01046 typedef FieldLayout<D> Layout_t;
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071 #ifdef PRINT_DEBUG
01072 msg << "Starting BC Calculation for face "
01073 << getFace() << "." << endl;
01074 #endif
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084 const NDIndex<D>& domain(A.getDomain());
01085
01086 NDIndex<D> dest_slab = AddGuardCells(domain,A.getGuardCellSizes());
01087
01088
01089
01090
01091 unsigned d = this->getFace()/2;
01092
01093 int offset;
01094
01095 CalcParallelPeriodicDomain(A,*this,dest_slab,offset);
01096
01097 Domain_t src_slab = dest_slab;
01098 src_slab[d] = src_slab[d] + offset;
01099
01100 #ifdef PRINT_DEBUG
01101 msg << "dest_slab = " << dest_slab << endl;
01102 msg << "src_slab = " << src_slab << endl;
01103
01104 #endif
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120 typedef multimap<Domain_t,LField_t*,less<Domain_t> > ReceiveMap_t;
01121
01122
01123
01124 TAU_PROFILE_START(findrectimer);
01125
01126 ReceiveMap_t receive_map;
01127
01128 TAU_PROFILE_STOP(findrectimer);
01129
01130
01131
01132 int receive_count = 0;
01133 int send_count = 0;
01134
01135
01136
01137 int bc_comm_tag;
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152 typedef vector<LField_t*> DestList_t;
01153 typedef vector<LField_t*> SrcList_t;
01154 typedef typename DestList_t::iterator DestListIterator_t;
01155 typedef typename SrcList_t::iterator SrcListIterator_t;
01156
01157 DestList_t dest_list;
01158 SrcList_t src_list;
01159
01160 dest_list.reserve(1);
01161 src_list.reserve(1);
01162
01163 typename Field_t::iterator_if lf_i;
01164
01165 #ifdef PRINT_DEBUG
01166 msg << "Starting dest & src domain calculation." << endl;
01167 #endif
01168
01169 for (lf_i = A.begin_if(); lf_i != A.end_if(); ++lf_i)
01170 {
01171 LField_t &lf = *lf_i->second;
01172
01173
01174
01175
01176 const Domain_t &lf_allocated = lf.getAllocated();
01177
01178 #ifdef PRINT_DEBUG
01179 msg << " Processing subdomain : " << lf_allocated << endl;
01180
01181 #endif
01182
01183 if (lf_allocated.touches(dest_slab))
01184 dest_list.push_back(&lf);
01185
01186
01187
01188
01189
01190 const Domain_t &lf_owned = lf.getOwned();
01191
01192 if (lf_owned.touches(src_slab))
01193 src_list.push_back(&lf);
01194 }
01195
01196 #ifdef PRINT_DEBUG
01197 msg << " dest_list has " << dest_list.size() << " elements." << endl;
01198 msg << " src_list has " << src_list.size() << " elements." << endl;
01199 #endif
01200
01201 DestListIterator_t dest_begin = dest_list.begin();
01202 DestListIterator_t dest_end = dest_list.end();
01203 SrcListIterator_t src_begin = src_list.begin();
01204 SrcListIterator_t src_end = src_list.end();
01205
01206
01207
01208 const Layout_t &layout = A.getLayout();
01209 const GuardCellSizes<D> &gc = A.getGuardCellSizes();
01210
01211 int nprocs = Ippl::getNodes();
01212
01213 if (nprocs > 1)
01214 {
01215 TAU_PROFILE_START(findrectimer);
01216
01217 #ifdef PRINT_DEBUG
01218 msg << "Starting receive calculation." << endl;
01219
01220 #endif
01221
01222
01223
01224
01225
01226
01227
01228 vector<bool> receive_mask(nprocs,false);
01229
01230 DestListIterator_t dest_i;
01231
01232 for (dest_i = dest_begin; dest_i != dest_end; ++dest_i)
01233 {
01234
01235
01236 LField_t &dest_lf = **dest_i;
01237
01238 const Domain_t &dest_lf_alloc = dest_lf.getAllocated();
01239
01240
01241
01242
01243
01244 const Domain_t dest_domain = dest_lf_alloc.intersect(dest_slab);
01245
01246 Domain_t src_domain = dest_domain;
01247 src_domain[d] = src_domain[d] + offset;
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260 typename Layout_t::touch_range_dv
01261 src_range(layout.touch_range_rdv(src_domain));
01262
01263
01264
01265
01266
01267
01268 typename Layout_t::touch_iterator_dv rv_i;
01269
01270 for (rv_i = src_range.first; rv_i != src_range.second; ++rv_i)
01271 {
01272
01273
01274
01275
01276
01277 const Domain_t rv_alloc = AddGuardCells(rv_i->first,gc);
01278
01279 Domain_t hit = src_domain.intersect(rv_alloc);
01280 hit[d] = hit[d] - offset;
01281
01282
01283
01284 typedef typename ReceiveMap_t::value_type value_type;
01285
01286 receive_map.insert(value_type(hit,&dest_lf));
01287
01288 #ifdef PRINT_DEBUG
01289 msg << " Need remote data for domain: " << hit << endl;
01290 #endif
01291
01292
01293
01294 int rnode = rv_i->second->getNode();
01295
01296 receive_mask[rnode] = true;
01297
01298 }
01299 }
01300
01301 receive_count = 0;
01302
01303 for (int iproc = 0; iproc < nprocs; ++iproc)
01304 if (receive_mask[iproc]) ++receive_count;
01305
01306
01307 #ifdef PRINT_DEBUG
01308 msg << " Expecting to see " << receive_count << " messages." << endl;
01309 msg << "Done with receive calculation." << endl;
01310
01311 #endif
01312
01313 TAU_PROFILE_STOP(findrectimer);
01314
01315
01316 TAU_PROFILE_START(sendtimer);
01317
01318 #ifdef PRINT_DEBUG
01319 msg << "Starting send calculation" << endl;
01320 #endif
01321
01322
01323
01324
01325
01326
01327
01328 vector<Message *> messages(nprocs);
01329 for (int miter=0; miter < nprocs; messages[miter++] = 0);
01330
01331
01332
01333 #ifdef PRINT_DEBUG
01334
01335
01336
01337 vector<int> ndomains(nprocs);
01338 for(int i = 0; i < nprocs; ++i) ndomains[i] = 0;
01339 #endif
01340
01341 SrcListIterator_t src_i;
01342
01343 for (src_i = src_begin; src_i != src_end; ++src_i)
01344 {
01345
01346
01347 LField_t &src_lf = **src_i;
01348
01349
01350
01351
01352
01353
01354
01355 const Domain_t &src_lf_owned = src_lf.getOwned();
01356 const Domain_t &src_lf_alloc = src_lf.getAllocated();
01357
01358
01359
01360
01361
01362 const Domain_t src_owned = src_lf_owned.intersect(src_slab);
01363 const Domain_t src_alloc = src_lf_alloc.intersect(src_slab);
01364
01365 Domain_t dest_owned = src_owned;
01366 dest_owned[d] = dest_owned[d] - offset;
01367
01368 Domain_t dest_alloc = src_alloc;
01369 dest_alloc[d] = dest_alloc[d] - offset;
01370
01371 #ifdef PRINT_DEBUG
01372 msg << " Considering LField with the domains:" << endl;
01373 msg << " owned = " << src_lf_owned << endl;
01374 msg << " alloc = " << src_lf_alloc << endl;
01375 msg << " The intersections with src_slab are:" << endl;
01376 msg << " owned = " << src_owned << endl;
01377 msg << " alloc = " << src_alloc << endl;
01378 #endif
01379
01380
01381
01382
01383 typename Layout_t::touch_range_dv
01384 dest_range(layout.touch_range_rdv(dest_owned,gc));
01385
01386 typename Layout_t::touch_iterator_dv rv_i;
01387
01388 #ifdef PRINT_DEBUG
01389 msg << " Touch calculation found "
01390 << distance(dest_range.first, dest_range.second)
01391 << " elements." << endl;
01392 #endif
01393
01394
01395 for (rv_i = dest_range.first; rv_i != dest_range.second; ++rv_i)
01396 {
01397
01398
01399
01400
01401 Domain_t hit = dest_alloc.intersect(rv_i->first);
01402 hit[d] = hit[d] + offset;
01403
01404
01405
01406 int rnode = rv_i->second->getNode();
01407
01408 #ifdef PRINT_DEBUG
01409 msg << " Overlap domain = " << rv_i->first << endl;
01410 msg << " Inters. domain = " << hit;
01411 msg << " --> node " << rnode << endl;
01412 #endif
01413
01414
01415
01416
01417
01418
01419
01420 Element_t compressed_value;
01421 LFI_t msgval = src_lf.begin(hit, compressed_value);
01422 msgval.TryCompress();
01423
01424
01425
01426 if (!messages[rnode])
01427 {
01428 messages[rnode] = new Message;
01429 PAssert(messages[rnode]);
01430 }
01431
01432 messages[rnode]->put(hit);
01433 messages[rnode]->put(msgval);
01434
01435 #ifdef PRINT_DEBUG
01436 ndomains[rnode]++;
01437 #endif
01438 }
01439 }
01440
01441
01442
01443 bc_comm_tag =
01444 Ippl::Comm->next_tag(BC_PARALLEL_PERIODIC_TAG,BC_TAG_CYCLE);
01445
01446 TAU_PROFILE_START(sendcommtimer);
01447
01448
01449
01450 for (int iproc = 0; iproc < nprocs; ++iproc)
01451 {
01452 if (messages[iproc])
01453 {
01454
01455 #ifdef PRINT_DEBUG
01456 msg << " ParallelPeriodicBCApply: Sending message to node "
01457 << iproc << endl
01458 << " number of domains = " << ndomains[iproc] << endl
01459 << " number of MsgItems = "
01460 << messages[iproc]->size() << endl;
01461 #endif
01462
01463 Ippl::Comm->send(messages[iproc], iproc, bc_comm_tag);
01464 ++send_count;
01465
01466 }
01467
01468 }
01469
01470 #ifdef PRINT_DEBUG
01471 msg << " Sent " << send_count << " messages" << endl;
01472 msg << "Done with send." << endl;
01473 #endif
01474
01475 TAU_PROFILE_STOP(sendcommtimer);
01476
01477 TAU_PROFILE_STOP(sendtimer);
01478
01479 }
01480
01481
01482 TAU_PROFILE_START(localstimer);
01483
01484
01485
01486
01487
01488 #ifdef PRINT_DEBUG
01489 msg << "Starting local calculation." << endl;
01490 #endif
01491
01492 DestListIterator_t dest_i;
01493
01494 for (dest_i = dest_begin; dest_i != dest_end; ++dest_i)
01495 {
01496
01497
01498 LField_t &dest_lf = **dest_i;
01499
01500 const Domain_t &dest_lf_alloc = dest_lf.getAllocated();
01501
01502 const Domain_t dest_domain = dest_lf_alloc.intersect(dest_slab);
01503
01504 Domain_t src_domain = dest_domain;
01505 src_domain[d] = src_domain[d] + offset;
01506
01507 SrcListIterator_t src_i;
01508
01509 for (src_i = src_begin; src_i != src_end; ++src_i)
01510 {
01511
01512
01513 LField_t &src_lf = **src_i;
01514
01515
01516
01517
01518
01519 const Domain_t &src_lf_owned = src_lf.getOwned();
01520 const Domain_t &src_lf_alloc = src_lf.getAllocated();
01521
01522
01523
01524
01525 if (src_domain.touches(src_lf_owned))
01526 {
01527
01528
01529
01530
01531 dest_lf.Uncompress();
01532
01533
01534
01535 Domain_t real_src_domain =
01536 src_domain.intersect(src_lf_alloc);
01537
01538 Domain_t real_dest_domain = real_src_domain;
01539 real_dest_domain[d] = real_dest_domain[d] - offset;
01540
01541 #ifdef PRINT_DEBUG
01542 msg << " Copying local data . . ." << endl;
01543 msg << " source domain = " << real_src_domain << endl;
01544 msg << " dest domain = " << real_dest_domain << endl;
01545 #endif
01546
01547
01548
01549 LFI_t lhs = dest_lf.begin(real_dest_domain);
01550 LFI_t rhs = src_lf.begin(real_src_domain);
01551
01552
01553
01554 if (this->getComponent() == BCondBase<T,D,M,C>::allComponents)
01555 {
01556 BrickExpression<D,LFI_t,LFI_t,OpPeriodic<T> >
01557 (lhs,rhs,OpPeriodic<T>()).apply();
01558 }
01559 else
01560 {
01561 BrickExpression<D,LFI_t,LFI_t,OpPeriodicComponent<T> >
01562 (lhs,rhs,OpPeriodicComponent<T>(this->getComponent())).apply();
01563 }
01564 }
01565 }
01566 }
01567
01568 #ifdef PRINT_DEBUG
01569 msg << "Done with local calculation." << endl;
01570 #endif
01571
01572 TAU_PROFILE_STOP(localstimer);
01573
01574
01575
01576
01577
01578
01579 if (nprocs > 1)
01580 {
01581
01582 TAU_PROFILE_START(rectimer);
01583
01584 #ifdef PRINT_DEBUG
01585 msg << "Starting receive..." << endl;
01586
01587 #endif
01588
01589 while (receive_count > 0)
01590 {
01591
01592
01593
01594 int any_node = COMM_ANY_NODE;
01595
01596 TAU_PROFILE_START(reccommtimer);
01597
01598 Message* message =
01599 Ippl::Comm->receive_block(any_node,bc_comm_tag);
01600 PAssert(message);
01601
01602 --receive_count;
01603
01604 TAU_PROFILE_STOP(reccommtimer);
01605
01606
01607
01608 int ndomains = message->size() / (D + 3);
01609
01610 #ifdef PRINT_DEBUG
01611 msg << " Message received from node "
01612 << any_node << "," << endl
01613 << " number of domains = " << ndomains << endl;
01614 #endif
01615
01616 for (int idomain=0; idomain < ndomains; ++idomain)
01617 {
01618
01619
01620
01621 Domain_t intersection;
01622 intersection.getMessage(*message);
01623 intersection[d] = intersection[d] - offset;
01624
01625
01626
01627 Element_t compressed_value;
01628 LFI_t rhs(compressed_value);
01629 rhs.getMessage(*message);
01630
01631 #ifdef PRINT_DEBUG
01632 msg << " Received remote overlap region = "
01633 << intersection << endl;
01634 #endif
01635
01636
01637
01638 typename ReceiveMap_t::iterator hit =
01639 receive_map.find(intersection);
01640
01641 PAssert(hit != receive_map.end());
01642
01643
01644
01645
01646
01647
01648
01649 LField<T,D> &lf = *(*hit).second;
01650
01651
01652
01653 #ifdef PRINT_DEBUG
01654 msg << " LHS compressed ? " << lf.IsCompressed();
01655 msg << ", LHS value = " << *lf.begin() << endl;
01656 msg << " RHS compressed ? " << rhs.IsCompressed();
01657 msg << ", RHS value = " << *rhs << endl;
01658 msg << " *rhs == *lf.begin() ? "
01659 << (*rhs == *lf.begin()) << endl;
01660 #endif
01661 if (!(rhs.IsCompressed() && lf.IsCompressed() &&
01662 (*rhs == *lf.begin())))
01663 {
01664
01665
01666 lf.Uncompress();
01667 LFI_t lhs = lf.begin(intersection);
01668
01669
01670
01671 #ifdef PRINT_DEBUG
01672 msg << " Doing copy." << endl;
01673 #endif
01674
01675 BrickExpression<D,LFI_t,LFI_t,OpAssign>(lhs,rhs).apply();
01676 }
01677
01678
01679
01680 receive_map.erase(hit);
01681 }
01682
01683 delete message;
01684 }
01685 TAU_PROFILE_STOP(rectimer);
01686
01687 #ifdef PRINT_DEBUG
01688 msg << "Done with receive." << endl;
01689 #endif
01690
01691 }
01692 }
01693
01694
01696
01697
01698
01699
01700 template<class T>
01701 struct OpExtrapolate
01702 {
01703 OpExtrapolate(const T& o, const T& s) : Offset(o), Slope(s) {}
01704 T Offset, Slope;
01705 };
01706 template<class T>
01707 inline void PETE_apply(const OpExtrapolate<T>& e, T& a, const T& b)
01708 { a = b*e.Slope + e.Offset; }
01709
01710
01711 template<class T>
01712 struct OpExtrapolateComponent
01713 {
01714 OpExtrapolateComponent(const T& o, const T& s, int c)
01715 : Offset(o), Slope(s), Component(c) {}
01716 T Offset, Slope;
01717 int Component;
01718 };
01719 template<class T>
01720 inline void PETE_apply(const OpExtrapolateComponent<T>& e, T& a, const T& b)
01721 {
01722 a[e.Component] = b[e.Component]*e.Slope[e.Component] + e.Offset[e.Component];
01723 }
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,char)
01741 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,bool)
01742 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,int)
01743 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,unsigned)
01744 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,short)
01745 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,long)
01746 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,float)
01747 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,double)
01748 COMPONENT_APPLY_BUILTIN(OpExtrapolateComponent,dcomplex)
01749
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01763
01764 template<class T, unsigned D, class M>
01765 void ExtrapolateFaceBCApply(ExtrapolateFace<T,D,M,Cell>& ef,
01766 Field<T,D,M,Cell>& A );
01767 template<class T, unsigned D, class M>
01768 void ExtrapolateFaceBCApply(ExtrapolateFace<T,D,M,Vert>& ef,
01769 Field<T,D,M,Vert>& A );
01770 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
01771 void ExtrapolateFaceBCApply(ExtrapolateFace<T,D,M,
01772 CartesianCentering<CE,D,NC> >& ef,
01773 Field<T,D,M,CartesianCentering<CE,D,NC> >& A );
01774
01775 template<class T, unsigned D, class M, class C>
01776 void ExtrapolateFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
01777 {
01778 ExtrapolateFaceBCApply(*this, A);
01779 }
01780
01781
01782 template<class T, unsigned D, class M, class C>
01783 void ExtrapolateFaceBCApply2(const NDIndex<D> &dest, const NDIndex<D> &src,
01784 LField<T,D> &fill, LField<T,D> &from, const NDIndex<D> &from_alloc,
01785 ExtrapolateFace<T,D,M,C> &ef)
01786 {
01787
01788
01789
01790
01791 if (fill.IsCompressed() && from.IsCompressed())
01792 {
01793
01794
01795
01796 T a = fill.getCompressedData(), aref = a;
01797 T b = from.getCompressedData();
01798 if (ef.getComponent() == BCondBase<T,D,M,C>::allComponents)
01799 {
01800 OpExtrapolate<T> op(ef.getOffset(),ef.getSlope());
01801 PETE_apply(op, a, b);
01802 }
01803 else
01804 {
01805 OpExtrapolateComponent<T>
01806 op(ef.getOffset(),ef.getSlope(),ef.getComponent());
01807 PETE_apply(op, a, b);
01808 }
01809 if (a == aref)
01810 {
01811
01812
01813 return;
01814 }
01815 }
01816
01817
01818
01819
01820 fill.Uncompress();
01821
01822 NDIndex<D> from_it = src.intersect(from_alloc);
01823 NDIndex<D> fill_it = dest.plugBase(from_it);
01824
01825
01826
01827 typedef typename LField<T,D>::iterator LFI;
01828 LFI lhs = fill.begin(fill_it);
01829 LFI rhs = from.begin(from_it);
01830
01831
01832
01833 if (ef.getComponent() == BCondBase<T,D,M,C>::allComponents)
01834 {
01835 BrickExpression<D,LFI,LFI,OpExtrapolate<T> >
01836 (lhs,rhs,OpExtrapolate<T>(ef.getOffset(),ef.getSlope())).apply();
01837 }
01838 else
01839 {
01840 BrickExpression<D,LFI,LFI,OpExtrapolateComponent<T> >
01841 (lhs,rhs,OpExtrapolateComponent<T>
01842 (ef.getOffset(),ef.getSlope(),ef.getComponent())).apply();
01843 }
01844 }
01845
01846
01847
01848
01849
01850
01851
01852 template<class T, unsigned D, class M>
01853 void ExtrapolateFaceBCApply(ExtrapolateFace<T,D,M,Cell>& ef,
01854 Field<T,D,M,Cell>& A )
01855 {
01856 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
01857 TAU_PROFILE("ExtrapolateFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
01858
01859
01860
01861
01862
01863 const NDIndex<D>& domain( A.getDomain() );
01864 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
01865
01866
01867
01868
01869
01870 unsigned d = ef.getFace()/2;
01871 int offset;
01872
01873
01874
01875
01876
01877
01878 if (ef.getFace() & 1)
01879 {
01880
01881
01882
01883
01884 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
01885
01886
01887
01888
01889 offset = 2*domain[d].max() + 1;
01890 }
01891 else
01892 {
01893
01894
01895
01896
01897 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
01898
01899
01900
01901
01902 offset = 2*domain[d].min() - 1;
01903 }
01904
01905
01906
01907 typename Field<T,D,M,Cell>::iterator_if fill_i;
01908 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
01909 {
01910
01911
01912
01913 LField<T,D> &fill = *(*fill_i).second;
01914
01915
01916
01917 const NDIndex<D> &fill_alloc = fill.getAllocated();
01918
01919
01920
01921
01922
01923
01924 if (slab.touches(fill_alloc))
01925 {
01926
01927
01928 NDIndex<D> dest = slab.intersect(fill_alloc);
01929
01930
01931
01932
01933
01934
01935
01936 NDIndex<D> src = dest;
01937
01938
01939
01940
01941 src[d] = offset - src[d];
01942
01943
01944
01945
01946 if (fill_alloc.contains(src))
01947 {
01948
01949
01950 ExtrapolateFaceBCApply2(dest, src, fill, fill,
01951 fill_alloc, ef);
01952 }
01953 else
01954 {
01955
01956
01957
01958 typename Field<T,D,M,Cell>::iterator_if from_i;
01959 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
01960 {
01961
01962
01963 LField<T,D> &from = *(*from_i).second;
01964 const NDIndex<D> &from_owned = from.getOwned();
01965 const NDIndex<D> &from_alloc = from.getAllocated();
01966
01967
01968
01969 if (src.touches(from_owned))
01970 ExtrapolateFaceBCApply2(dest, src, fill, from,
01971 from_alloc, ef);
01972 }
01973 }
01974 }
01975 }
01976 }
01977
01978
01979
01980
01981
01982
01983
01984 template<class T, unsigned D, class M>
01985 void ExtrapolateFaceBCApply(ExtrapolateFace<T,D,M,Vert>& ef,
01986 Field<T,D,M,Vert>& A )
01987 {
01988 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
01989 TAU_PROFILE("ExtrapolateFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
01990
01991
01992
01993
01994
01995 const NDIndex<D>& domain(A.getDomain());
01996 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
01997
01998
01999
02000
02001
02002 unsigned d = ef.getFace()/2;
02003 int offset;
02004
02005
02006
02007
02008
02009
02010 if ( ef.getFace() & 1 )
02011 {
02012
02013
02014
02015
02016 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
02017
02018
02019
02020
02021
02022
02023 offset = 2*domain[d].max() + 1 - 1;
02024 }
02025 else
02026 {
02027
02028
02029
02030
02031 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
02032
02033
02034
02035
02036
02037 offset = 2*domain[d].min() - 1 + 1;
02038 }
02039
02040
02041
02042 typename Field<T,D,M,Vert>::iterator_if fill_i;
02043 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
02044 {
02045
02046
02047
02048 LField<T,D> &fill = *(*fill_i).second;
02049
02050
02051 const NDIndex<D> &fill_alloc = fill.getAllocated();
02052
02053
02054
02055
02056
02057 if ( slab.touches( fill_alloc ) )
02058 {
02059
02060
02061 NDIndex<D> dest = slab.intersect( fill_alloc );
02062
02063
02064
02065
02066
02067
02068
02069 NDIndex<D> src = dest;
02070
02071
02072
02073
02074 src[d] = offset - src[d];
02075
02076
02077
02078
02079 if (fill_alloc.contains(src))
02080 {
02081
02082
02083 ExtrapolateFaceBCApply2(dest, src, fill, fill,
02084 fill_alloc, ef);
02085 }
02086 else
02087 {
02088
02089
02090
02091 typename Field<T,D,M,Vert>::iterator_if from_i;
02092 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
02093 {
02094
02095
02096 LField<T,D> &from = *(*from_i).second;
02097 const NDIndex<D> &from_owned = from.getOwned();
02098 const NDIndex<D> &from_alloc = from.getAllocated();
02099
02100
02101
02102 if (src.touches(from_owned))
02103 ExtrapolateFaceBCApply2(dest, src, fill, from,
02104 from_alloc, ef);
02105 }
02106 }
02107 }
02108 }
02109 }
02110
02111
02112
02113
02114
02115
02116 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
02117 void ExtrapolateFaceBCApply(ExtrapolateFace<T,D,M,
02118 CartesianCentering<CE,D,NC> >& ef,
02119 Field<T,D,M,CartesianCentering<CE,D,NC> >& A )
02120 {
02121 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
02122 TAU_PROFILE("ExtrapolateFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
02123
02124
02125
02126
02127
02128 const NDIndex<D>& domain( A.getDomain() );
02129 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
02130
02131
02132
02133
02134
02135 unsigned d = ef.getFace()/2;
02136 int offset;
02137
02138
02139
02140
02141
02142
02143 if ( ef.getFace() & 1 )
02144 {
02145
02146
02147
02148
02149
02150 if (ef.getComponent() == BCondBase<T,D,M,CartesianCentering<CE,D,NC> >::
02151 allComponents)
02152 {
02153
02154
02155 CenteringEnum centering0 = CE[0 + d*NC];
02156 for (int c=1; c<NC; c++)
02157 {
02158
02159 if (CE[c + d*NC] != centering0)
02160 ERRORMSG("ExtrapolateFaceBCApply: BCond thinks all components"
02161 << " have same centering along direction " << d
02162 << ", but it isn't so." << endl);
02163 }
02164
02165
02166
02167
02168
02169
02170
02171 slab[d] = Index(domain[d].max() + 1,
02172 domain[d].max() + A.rightGuard(d));
02173
02174 if (centering0 == CELL)
02175 {
02176 offset = 2*domain[d].max() + 1 ;
02177 }
02178 else
02179 {
02180 offset = 2*domain[d].max() + 1 - 1;
02181 }
02182 }
02183 else
02184 {
02185
02186
02187 if (CE[ef.getComponent() + d*NC] == CELL)
02188 {
02189
02190
02191
02192 int highcell = A.get_mesh().gridSizes[d] - 2;
02193 slab[d] = Index(highcell + 1, highcell + A.rightGuard(d));
02194
02195
02196 offset = 2*highcell + 1 ;
02197 }
02198 else
02199 {
02200
02201
02202
02203 int highvert = A.get_mesh().gridSizes[d] - 1;
02204 slab[d] = Index(highvert + 1, highvert + A.rightGuard(d));
02205
02206
02207 offset = 2*highvert + 1 - 1;
02208 }
02209 }
02210 }
02211 else
02212 {
02213
02214
02215
02216
02217 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
02218
02219
02220
02221
02222
02223
02224 if (ef.getComponent() == BCondBase<T,D,M,CartesianCentering<CE,D,NC> >::
02225 allComponents)
02226 {
02227
02228
02229 CenteringEnum centering0 = CE[0 + d*NC];
02230 for (int c=1; c<NC; c++)
02231 {
02232
02233
02234 if (CE[c + d*NC] != centering0)
02235 ERRORMSG("ExtrapolateFaceBCApply: BCond thinks all components"
02236 << " have same centering along direction " << d
02237 << ", but it isn't so." << endl);
02238 }
02239
02240
02241
02242
02243 if (centering0 == CELL)
02244 {
02245 offset = 2*domain[d].min() - 1;
02246 }
02247 else
02248 {
02249 offset = 2*domain[d].min() - 1 + 1;
02250 }
02251 }
02252 else
02253 {
02254
02255
02256
02257 if (CE[ef.getComponent() + d*NC] == CELL)
02258 {
02259 offset = 2*domain[d].min() - 1;
02260 }
02261 else
02262 {
02263 offset = 2*domain[d].min() - 1 + 1;
02264 }
02265 }
02266 }
02267
02268
02269
02270 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if fill_i;
02271 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
02272 {
02273
02274
02275
02276 LField<T,D> &fill = *(*fill_i).second;
02277
02278
02279
02280 const NDIndex<D> &fill_alloc = fill.getAllocated();
02281
02282
02283
02284
02285
02286
02287 if ( slab.touches( fill_alloc ) )
02288 {
02289
02290
02291 NDIndex<D> dest = slab.intersect( fill_alloc );
02292
02293
02294
02295
02296
02297
02298
02299 NDIndex<D> src = dest;
02300
02301
02302
02303
02304 src[d] = offset - src[d];
02305
02306
02307
02308
02309 if (fill_alloc.contains(src))
02310 {
02311
02312
02313 ExtrapolateFaceBCApply2(dest, src, fill, fill,
02314 fill_alloc, ef);
02315 }
02316 else
02317 {
02318
02319
02320
02321 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if
02322 from_i;
02323 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
02324 {
02325
02326
02327 LField<T,D> &from = *(*from_i).second;
02328 const NDIndex<D> &from_owned = from.getOwned();
02329 const NDIndex<D> &from_alloc = from.getAllocated();
02330
02331
02332
02333 if (src.touches(from_owned))
02334 ExtrapolateFaceBCApply2(dest, src, fill, from,
02335 from_alloc, ef);
02336 }
02337 }
02338 }
02339 }
02340 }
02341
02342
02343
02345
02346
02347
02348
02349 template<class T>
02350 struct OpExtrapolateAndZero
02351 {
02352 OpExtrapolateAndZero(const T& o, const T& s) : Offset(o), Slope(s) {}
02353 T Offset, Slope;
02354 };
02355 template<class T>
02356 inline void PETE_apply(const OpExtrapolateAndZero<T>& e, T& a, const T& b)
02357 { a = b*e.Slope + e.Offset; }
02358
02359
02360 template<class T>
02361 struct OpExtrapolateAndZeroComponent
02362 {
02363 OpExtrapolateAndZeroComponent(const T& o, const T& s, int c)
02364 : Offset(o), Slope(s), Component(c) {}
02365 T Offset, Slope;
02366 int Component;
02367 };
02368 template<class T>
02369 inline void PETE_apply(const OpExtrapolateAndZeroComponent<T>& e, T& a,
02370 const T& b)
02371 {
02372 a[e.Component] = b[e.Component]*e.Slope[e.Component] + e.Offset[e.Component];
02373 }
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,char)
02393 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,bool)
02394 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,int)
02395 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,unsigned)
02396 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,short)
02397 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,long)
02398 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,float)
02399 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,double)
02400 COMPONENT_APPLY_BUILTIN(OpExtrapolateAndZeroComponent,dcomplex)
02401
02402
02403 template<class T>
02404 struct OpAssignComponent
02405 {
02406 OpAssignComponent(int c)
02407 : Component(c) { }
02408 int Component;
02409 };
02410
02411 template<class T, class T1>
02412 inline void PETE_apply(const OpAssignComponent<T>& e, T& a, const T1& b)
02413 {
02414 a[e.Component] = b;
02415 }
02416
02417 COMPONENT_APPLY_BUILTIN(OpAssignComponent,char)
02418 COMPONENT_APPLY_BUILTIN(OpAssignComponent,bool)
02419 COMPONENT_APPLY_BUILTIN(OpAssignComponent,int)
02420 COMPONENT_APPLY_BUILTIN(OpAssignComponent,unsigned)
02421 COMPONENT_APPLY_BUILTIN(OpAssignComponent,short)
02422 COMPONENT_APPLY_BUILTIN(OpAssignComponent,long)
02423 COMPONENT_APPLY_BUILTIN(OpAssignComponent,float)
02424 COMPONENT_APPLY_BUILTIN(OpAssignComponent,double)
02425 COMPONENT_APPLY_BUILTIN(OpAssignComponent,dcomplex)
02426
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02440
02441 template<class T, unsigned D, class M>
02442 void ExtrapolateAndZeroFaceBCApply(ExtrapolateAndZeroFace<T,D,M,Cell>& ef,
02443 Field<T,D,M,Cell>& A );
02444 template<class T, unsigned D, class M>
02445 void ExtrapolateAndZeroFaceBCApply(ExtrapolateAndZeroFace<T,D,M,Vert>& ef,
02446 Field<T,D,M,Vert>& A );
02447 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
02448 void ExtrapolateAndZeroFaceBCApply(ExtrapolateAndZeroFace<T,D,M,
02449 CartesianCentering<CE,D,NC> >& ef,
02450 Field<T,D,M,CartesianCentering<CE,D,NC> >& A );
02451
02452 template<class T, unsigned D, class M, class C>
02453 void ExtrapolateAndZeroFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
02454 {
02455 ExtrapolateAndZeroFaceBCApply(*this, A);
02456 }
02457
02458
02459 template<class T, unsigned D, class M, class C>
02460 void ExtrapolateAndZeroFaceBCApply2(const NDIndex<D> &dest,
02461 const NDIndex<D> &src, LField<T,D> &fill, LField<T,D> &from,
02462 const NDIndex<D> &from_alloc, ExtrapolateAndZeroFace<T,D,M,C> &ef)
02463 {
02464
02465
02466
02467
02468 if (fill.IsCompressed() && from.IsCompressed())
02469 {
02470
02471
02472
02473 T a = fill.getCompressedData(), aref = a;
02474 T b = from.getCompressedData();
02475 if (ef.getComponent() == BCondBase<T,D,M,C>::allComponents)
02476 {
02477 OpExtrapolateAndZero<T> op(ef.getOffset(),ef.getSlope());
02478 PETE_apply(op, a, b);
02479 }
02480 else
02481 {
02482 OpExtrapolateAndZeroComponent<T>
02483 op(ef.getOffset(),ef.getSlope(),ef.getComponent());
02484 PETE_apply(op, a, b);
02485 }
02486 if (a == aref)
02487 {
02488
02489
02490 return;
02491 }
02492 }
02493
02494
02495
02496
02497 fill.Uncompress();
02498
02499 NDIndex<D> from_it = src.intersect(from_alloc);
02500 NDIndex<D> fill_it = dest.plugBase(from_it);
02501
02502
02503
02504 typedef typename LField<T,D>::iterator LFI;
02505 LFI lhs = fill.begin(fill_it);
02506 LFI rhs = from.begin(from_it);
02507
02508
02509
02510 if (ef.getComponent() == BCondBase<T,D,M,C>::allComponents)
02511 {
02512 BrickExpression< D, LFI, LFI, OpExtrapolateAndZero<T> >
02513 (lhs, rhs,
02514 OpExtrapolateAndZero<T>(ef.getOffset(),ef.getSlope())).apply();
02515 }
02516 else
02517 {
02518 BrickExpression< D, LFI, LFI, OpExtrapolateAndZeroComponent<T> >
02519 (lhs, rhs,
02520 OpExtrapolateAndZeroComponent<T>
02521 (ef.getOffset(),ef.getSlope(),ef.getComponent())).apply();
02522 }
02523 }
02524
02525
02526 template<class T, unsigned D, class M, class C>
02527 void ExtrapolateAndZeroFaceBCApply3(const NDIndex<D> &dest,
02528 LField<T,D> &fill, ExtrapolateAndZeroFace<T,D,M,C> &ef)
02529 {
02530
02531
02532
02533
02534 if (fill.IsCompressed())
02535 {
02536
02537
02538
02539 if (ef.getComponent() == BCondBase<T,D,M,C>::allComponents)
02540 {
02541 if (fill.getCompressedData() == T(0))
02542 return;
02543 }
02544 else
02545 {
02546 typedef typename AppTypeTraits<T>::Element_t T1;
02547
02548
02549
02550
02551
02552 T a = fill.getCompressedData(), aref = a;
02553 OpAssignComponent<T> op(ef.getComponent());
02554 PETE_apply(op, a, T1(0));
02555 if (a == aref)
02556 return;
02557 }
02558 }
02559
02560
02561
02562
02563 fill.Uncompress();
02564
02565
02566
02567 typedef typename LField<T,D>::iterator LFI;
02568 LFI lhs = fill.begin(dest);
02569
02570
02571
02572 if (ef.getComponent() == BCondBase<T,D,M,C>::allComponents)
02573 {
02574 BrickExpression<D,LFI,PETE_Scalar<T>,OpAssign >
02575 (lhs,PETE_Scalar<T>(T(0)),OpAssign()).apply();
02576 }
02577 else
02578 {
02579 typedef typename AppTypeTraits<T>::Element_t T1;
02580
02581 BrickExpression<D,LFI,PETE_Scalar<T1>,OpAssignComponent<T> >
02582 (lhs,PETE_Scalar<T1>(T1(0)),OpAssignComponent<T>
02583 (ef.getComponent())).apply();
02584 }
02585 }
02586
02587
02588
02589
02590
02591
02592
02593
02594 template<class T, unsigned D, class M>
02595 void ExtrapolateAndZeroFaceBCApply(ExtrapolateAndZeroFace<T,D,M,Cell>& ef,
02596 Field<T,D,M,Cell>& A )
02597 {
02598 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
02599 TAU_PROFILE("ExtrapolateAndZeroFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
02600
02601
02602
02603
02604
02605 const NDIndex<D>& domain( A.getDomain() );
02606 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
02607
02608
02609
02610
02611
02612 unsigned d = ef.getFace()/2;
02613 int offset;
02614
02615
02616
02617
02618
02619
02620 if (ef.getFace() & 1)
02621 {
02622
02623
02624
02625
02626 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
02627
02628
02629
02630
02631 offset = 2*domain[d].max() + 1;
02632 }
02633 else
02634 {
02635
02636
02637
02638
02639 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
02640
02641
02642
02643
02644 offset = 2*domain[d].min() - 1;
02645 }
02646
02647
02648
02649 typename Field<T,D,M,Cell>::iterator_if fill_i;
02650 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
02651 {
02652
02653
02654
02655 LField<T,D> &fill = *(*fill_i).second;
02656
02657
02658
02659 const NDIndex<D> &fill_alloc = fill.getAllocated();
02660
02661
02662
02663
02664
02665
02666 if (slab.touches(fill_alloc))
02667 {
02668
02669
02670 NDIndex<D> dest = slab.intersect(fill_alloc);
02671
02672
02673
02674
02675
02676
02677
02678 NDIndex<D> src = dest;
02679
02680
02681
02682
02683 src[d] = offset - src[d];
02684
02685
02686
02687
02688 if (fill_alloc.contains(src))
02689 {
02690
02691
02692 ExtrapolateAndZeroFaceBCApply2(dest, src, fill, fill,
02693 fill_alloc, ef);
02694 }
02695 else
02696 {
02697
02698
02699
02700 typename Field<T,D,M,Cell>::iterator_if from_i;
02701 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
02702 {
02703
02704
02705 LField<T,D> &from = *(*from_i).second;
02706 const NDIndex<D> &from_owned = from.getOwned();
02707 const NDIndex<D> &from_alloc = from.getAllocated();
02708
02709
02710
02711 if (src.touches(from_owned))
02712 ExtrapolateAndZeroFaceBCApply2(dest, src, fill, from,
02713 from_alloc, ef);
02714 }
02715 }
02716 }
02717 }
02718 }
02719
02720
02721
02722
02723
02724
02725
02726
02727 template<class T, unsigned D, class M>
02728 void ExtrapolateAndZeroFaceBCApply(ExtrapolateAndZeroFace<T,D,M,Vert>& ef,
02729 Field<T,D,M,Vert>& A )
02730 {
02731 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
02732 TAU_PROFILE("ExtrapolateAndZeroFaceBCApply()", taustr, TAU_FIELD |
02733 TAU_ASSIGN);
02734
02735
02736
02737
02738
02739 const NDIndex<D>& domain(A.getDomain());
02740 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
02741
02742 NDIndex<D> phys = slab;
02743
02744
02745
02746
02747
02748 unsigned d = ef.getFace()/2;
02749 int offset;
02750
02751
02752
02753
02754
02755
02756 if ( ef.getFace() & 1 )
02757 {
02758
02759
02760
02761
02762 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
02763
02764
02765
02766 phys[d] = Index( domain[d].max(), domain[d].max(), 1);
02767
02768
02769
02770
02771
02772
02773 offset = 2*domain[d].max() + 1 - 1;
02774 }
02775 else
02776 {
02777
02778
02779
02780
02781 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
02782
02783
02784
02785 phys[d] = Index( domain[d].min(), domain[d].min(), 1);
02786
02787
02788
02789
02790
02791
02792 offset = 2*domain[d].min() - 1 + 1;
02793 }
02794
02795
02796
02797 typename Field<T,D,M,Vert>::iterator_if fill_i;
02798 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
02799 {
02800
02801
02802
02803 LField<T,D> &fill = *(*fill_i).second;
02804
02805
02806
02807
02808
02809 const NDIndex<D> &fill_alloc = fill.getAllocated();
02810
02811
02812 if (phys.touches(fill_alloc))
02813 {
02814
02815
02816
02817 NDIndex<D> dest = phys.intersect(fill_alloc);
02818
02819
02820
02821 ExtrapolateAndZeroFaceBCApply3(dest, fill, ef);
02822 }
02823
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833 if ( slab.touches( fill_alloc ) )
02834 {
02835
02836
02837 NDIndex<D> dest = slab.intersect( fill_alloc );
02838
02839
02840
02841
02842
02843
02844
02845 NDIndex<D> src = dest;
02846
02847
02848
02849
02850 src[d] = offset - src[d];
02851
02852
02853
02854
02855 if (fill_alloc.contains(src))
02856 {
02857
02858
02859 ExtrapolateAndZeroFaceBCApply2(dest, src, fill, fill,
02860 fill_alloc, ef);
02861 }
02862 else
02863 {
02864
02865
02866
02867 typename Field<T,D,M,Vert>::iterator_if from_i;
02868 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
02869 {
02870
02871
02872 LField<T,D> &from = *(*from_i).second;
02873 const NDIndex<D> &from_owned = from.getOwned();
02874 const NDIndex<D> &from_alloc = from.getAllocated();
02875
02876
02877
02878 if (src.touches(from_owned))
02879 ExtrapolateAndZeroFaceBCApply2(dest, src, fill, from,
02880 from_alloc, ef);
02881 }
02882 }
02883 }
02884 }
02885 }
02886
02887
02888
02889
02890
02891
02892
02893 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
02894 void ExtrapolateAndZeroFaceBCApply(ExtrapolateAndZeroFace<T,D,M,
02895 CartesianCentering<CE,D,NC> >& ef,
02896 Field<T,D,M,CartesianCentering<CE,D,NC> >& A )
02897 {
02898 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
02899 TAU_PROFILE("ExtrapolateAndZeroFaceBCApply()", taustr, TAU_FIELD |
02900 TAU_ASSIGN);
02901
02902
02903
02904
02905
02906 const NDIndex<D>& domain( A.getDomain() );
02907 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
02908
02909 NDIndex<D> phys = slab;
02910
02911
02912
02913
02914
02915 unsigned d = ef.getFace()/2;
02916 int offset;
02917 bool setPhys = false;
02918
02919
02920
02921
02922
02923
02924 if ( ef.getFace() & 1 )
02925 {
02926
02927
02928
02929
02930
02931 if (ef.getComponent() == BCondBase<T,D,M,CartesianCentering<CE,D,NC> >::
02932 allComponents)
02933 {
02934
02935
02936 CenteringEnum centering0 = CE[0 + d*NC];
02937 for (int c=1; c<NC; c++)
02938 {
02939
02940 if (CE[c + d*NC] != centering0)
02941 ERRORMSG(
02942 "ExtrapolateAndZeroFaceBCApply: BCond thinks all components"
02943 << " have same centering along direction " << d
02944 << ", but it isn't so." << endl);
02945 }
02946
02947
02948
02949
02950
02951
02952
02953 slab[d] = Index(domain[d].max() + 1,
02954 domain[d].max() + A.rightGuard(d));
02955
02956 if (centering0 == CELL)
02957 {
02958 offset = 2*domain[d].max() + 1 ;
02959 }
02960 else
02961 {
02962 offset = 2*domain[d].max() + 1 - 1;
02963
02964
02965
02966 phys[d] = Index( domain[d].max(), domain[d].max(), 1);
02967 setPhys = true;
02968 }
02969 }
02970 else
02971 {
02972
02973
02974 if (CE[ef.getComponent() + d*NC] == CELL)
02975 {
02976
02977
02978
02979 int highcell = A.get_mesh().gridSizes[d] - 2;
02980 slab[d] = Index(highcell + 1, highcell + A.rightGuard(d));
02981
02982
02983 offset = 2*highcell + 1 ;
02984 }
02985 else
02986 {
02987
02988
02989
02990
02991 int highvert = A.get_mesh().gridSizes[d] - 1;
02992 slab[d] = Index(highvert + 1, highvert + A.rightGuard(d));
02993
02994
02995
02996 offset = 2*highvert + 1 - 1;
02997
02998
02999
03000 phys[d] = Index( highvert, highvert, 1 );
03001 setPhys = true;
03002 }
03003 }
03004 }
03005 else
03006 {
03007
03008
03009
03010
03011 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03012
03013
03014
03015
03016
03017
03018 if (ef.getComponent() == BCondBase<T,D,M,CartesianCentering<CE,D,NC> >::
03019 allComponents)
03020 {
03021
03022
03023 CenteringEnum centering0 = CE[0 + d*NC];
03024 for (int c=1; c<NC; c++)
03025 {
03026
03027
03028 if (CE[c + d*NC] != centering0)
03029 ERRORMSG(
03030 "ExtrapolateAndZeroFaceBCApply: BCond thinks all components"
03031 << " have same centering along direction " << d
03032 << ", but it isn't so." << endl);
03033 }
03034
03035
03036
03037
03038 if (centering0 == CELL)
03039 {
03040 offset = 2*domain[d].min() - 1;
03041 }
03042 else
03043 {
03044 offset = 2*domain[d].min() - 1 + 1;
03045
03046
03047
03048 phys[d] = Index(domain[d].min(), domain[d].min(), 1);
03049 setPhys = true;
03050 }
03051 }
03052 else
03053 {
03054
03055
03056
03057 if (CE[ef.getComponent() + d*NC] == CELL)
03058 {
03059 offset = 2*domain[d].min() - 1;
03060 }
03061 else
03062 {
03063 offset = 2*domain[d].min() - 1 + 1;
03064
03065
03066
03067 phys[d] = Index(domain[d].min(), domain[d].min(), 1);
03068 setPhys = true;
03069 }
03070 }
03071 }
03072
03073
03074
03075 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if fill_i;
03076 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03077 {
03078
03079
03080
03081 LField<T,D> &fill = *(*fill_i).second;
03082
03083
03084
03085
03086 const NDIndex<D> &fill_alloc = fill.getAllocated();
03087
03088 if (setPhys)
03089 {
03090
03091
03092
03093 if (phys.touches(fill_alloc))
03094 {
03095
03096
03097
03098 NDIndex<D> dest = phys.intersect(fill_alloc);
03099
03100
03101
03102 ExtrapolateAndZeroFaceBCApply3(dest, fill, ef);
03103 }
03104 }
03105
03106
03107
03108
03109
03110
03111
03112
03113
03114
03115 if ( slab.touches( fill_alloc ) )
03116 {
03117
03118
03119 NDIndex<D> dest = slab.intersect( fill_alloc );
03120
03121
03122
03123
03124
03125
03126
03127 NDIndex<D> src = dest;
03128
03129
03130
03131
03132 src[d] = offset - src[d];
03133
03134
03135
03136
03137 if (fill_alloc.contains(src))
03138 {
03139
03140
03141 ExtrapolateAndZeroFaceBCApply2(dest, src, fill, fill,
03142 fill_alloc, ef);
03143 }
03144 else
03145 {
03146
03147
03148
03149 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if
03150 from_i;
03151 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03152 {
03153
03154
03155 LField<T,D> &from = *(*from_i).second;
03156 const NDIndex<D> &from_owned = from.getOwned();
03157 const NDIndex<D> &from_alloc = from.getAllocated();
03158
03159
03160
03161 if (src.touches(from_owned))
03162 ExtrapolateAndZeroFaceBCApply2(dest, src, fill, from,
03163 from_alloc, ef);
03164 }
03165 }
03166 }
03167 }
03168
03169 }
03170
03171
03172
03173
03175
03176
03177
03178
03179 template<class T>
03180 struct OpBCFunctionEq
03181 {
03182 OpBCFunctionEq(T (*func)(const T&)) : Func(func) {}
03183 T (*Func)(const T&);
03184 };
03185 template<class T>
03186
03187 inline void PETE_apply(const OpBCFunctionEq<T>& e, T& a, T& b)
03188 {
03189 a = e.Func(b);
03190 }
03191
03192
03193
03194
03196
03197
03198
03199
03200
03201
03202
03203
03204
03205
03206
03208
03209 template<class T, unsigned D, class M>
03210 void FunctionFaceBCApply(FunctionFace<T,D,M,Cell>& ff,
03211 Field<T,D,M,Cell>& A );
03212 template<class T, unsigned D, class M>
03213 void FunctionFaceBCApply(FunctionFace<T,D,M,Vert>& ff,
03214 Field<T,D,M,Vert>& A );
03215 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
03216 void FunctionFaceBCApply(FunctionFace<T,D,M,
03217 CartesianCentering<CE,D,NC> >& ff,
03218 Field<T,D,M,CartesianCentering<CE,D,NC> >& A );
03219
03220 template<class T, unsigned D, class M, class C>
03221 void FunctionFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
03222 {
03223 TAU_TYPE_STRING(taustr, "void (" + CT(A) + " )" );
03224 TAU_PROFILE("FunctionFace::apply()", taustr, TAU_FIELD | TAU_ASSIGN);
03225 FunctionFaceBCApply(*this, A);
03226 }
03227
03228
03229
03230
03231
03232 template<class T, unsigned D, class M>
03233 void FunctionFaceBCApply(FunctionFace<T,D,M,Cell>& ff,
03234 Field<T,D,M,Cell>& A )
03235 {
03236 TAU_TYPE_STRING(taustr, "void (" + CT(ff) + ", " + CT(A) + " )" );
03237 TAU_PROFILE("FunctionFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
03238
03239
03240
03241
03242
03243 const NDIndex<D>& domain( A.getDomain() );
03244 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
03245 unsigned d = ff.getFace()/2;
03246 int offset;
03247 if ( ff.getFace() & 1 ) {
03248
03249 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
03250 offset = 2*domain[d].max() + 1;
03251 }
03252 else {
03253 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03254 offset = 2*domain[d].min() - 1;
03255 }
03256
03257
03258 typename Field<T,D,M,Cell>::iterator_if fill_i;
03259 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03260 {
03261
03262 LField<T,D> &fill = *(*fill_i).second;
03263 const NDIndex<D> &fill_alloc = fill.getAllocated();
03264 if ( slab.touches( fill_alloc ) )
03265 {
03266
03267 NDIndex<D> dest = slab.intersect( fill_alloc );
03268
03269
03270 NDIndex<D> src = dest;
03271 src[d] = offset - src[d];
03272
03273
03274 typename Field<T,D,M,Cell>::iterator_if from_i;
03275 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03276 {
03277
03278 LField<T,D> &from = *(*from_i).second;
03279 const NDIndex<D> &from_owned = from.getOwned();
03280 const NDIndex<D> &from_alloc = from.getAllocated();
03281
03282 if ( src.touches( from_owned ) )
03283 {
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293
03294
03295
03296 fill.Uncompress();
03297
03298 NDIndex<D> from_it = src.intersect( from_alloc );
03299 NDIndex<D> fill_it = dest.plugBase( from_it );
03300
03301 typedef typename LField<T,D>::iterator LFI;
03302 LFI lhs = fill.begin(fill_it);
03303 LFI rhs = from.begin(from_it);
03304
03305 BrickExpression<D,LFI,LFI,OpBCFunctionEq<T> >
03306 (lhs,rhs,OpBCFunctionEq<T>(ff.Func)).apply();
03307 }
03308 }
03309 }
03310 }
03311 }
03312
03313
03314
03315
03316
03317 template<class T, unsigned D, class M>
03318 void FunctionFaceBCApply(FunctionFace<T,D,M,Vert>& ff,
03319 Field<T,D,M,Vert>& A )
03320 {
03321 TAU_TYPE_STRING(taustr, "void (" + CT(ff) + ", " + CT(A) + " )" );
03322 TAU_PROFILE("FunctionFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
03323
03324
03325
03326
03327
03328 const NDIndex<D>& domain( A.getDomain() );
03329 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
03330 unsigned d = ff.getFace()/2;
03331 int offset;
03332 if ( ff.getFace() & 1 ) {
03333
03334 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
03335
03336
03337 offset = 2*domain[d].max() + 1 - 1;
03338 }
03339 else {
03340 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03341
03342
03343 offset = 2*domain[d].min() - 1 + 1;
03344 }
03345
03346
03347 typename Field<T,D,M,Vert>::iterator_if fill_i;
03348 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03349 {
03350
03351 LField<T,D> &fill = *(*fill_i).second;
03352 const NDIndex<D> &fill_alloc = fill.getAllocated();
03353 if ( slab.touches( fill_alloc ) )
03354 {
03355
03356 NDIndex<D> dest = slab.intersect( fill_alloc );
03357
03358
03359 NDIndex<D> src = dest;
03360 src[d] = offset - src[d];
03361
03362
03363 typename Field<T,D,M,Vert>::iterator_if from_i;
03364 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03365 {
03366
03367 LField<T,D> &from = *(*from_i).second;
03368 const NDIndex<D> &from_owned = from.getOwned();
03369 const NDIndex<D> &from_alloc = from.getAllocated();
03370
03371 if ( src.touches( from_owned ) )
03372 {
03373
03374
03375
03376
03377
03378
03379
03380
03381
03382
03383
03384
03385 fill.Uncompress();
03386
03387 NDIndex<D> from_it = src.intersect( from_alloc );
03388 NDIndex<D> fill_it = dest.plugBase( from_it );
03389
03390 typedef typename LField<T,D>::iterator LFI;
03391 LFI lhs = fill.begin(fill_it);
03392 LFI rhs = from.begin(from_it);
03393
03394 BrickExpression<D,LFI,LFI,OpBCFunctionEq<T> >
03395 (lhs,rhs,OpBCFunctionEq<T>(ff.Func)).apply();
03396 }
03397 }
03398 }
03399 }
03400 }
03401
03402
03403
03404
03405
03406 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
03407 void FunctionFaceBCApply(FunctionFace<T,D,M,
03408 CartesianCentering<CE,D,NC> >& ff,
03409 Field<T,D,M,CartesianCentering<CE,D,NC> >& A )
03410 {
03411 TAU_TYPE_STRING(taustr, "void (" + CT(ff) + ", " + CT(A) + " )" );
03412 TAU_PROFILE("FunctionFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
03413
03414
03415
03416
03417
03418 const NDIndex<D>& domain( A.getDomain() );
03419 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
03420 unsigned d = ff.getFace()/2;
03421 int offset;
03422 if ( ff.getFace() & 1 ) {
03423
03424 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
03425
03426
03427 CenteringEnum centering0 = CE[0 + d*NC];
03428 for (int c=1; c<NC; c++) {
03429 if (CE[c + d*NC] != centering0)
03430 ERRORMSG("FunctionFaceBCApply: BCond thinks all components have"
03431 << " same centering along direction " << d
03432 << ", but it isn't so." << endl);
03433 }
03434
03435 if (centering0 == CELL) {
03436 offset = 2*domain[d].max() + 1;
03437 } else {
03438 offset = 2*domain[d].max() + 1 - 1;
03439 }
03440 }
03441 else {
03442 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03443
03444
03445 CenteringEnum centering0 = CE[0 + d*NC];
03446 for (int c=1; c<NC; c++) {
03447 if (CE[c + d*NC] != centering0)
03448 ERRORMSG("FunctionFaceBCApply: BCond thinks all components have"
03449 << " same centering along direction " << d
03450 << ", but it isn't so." << endl);
03451 }
03452
03453 if (centering0 == CELL) {
03454 offset = 2*domain[d].min() - 1;
03455 } else {
03456 offset = 2*domain[d].min() - 1 + 1;
03457 }
03458 }
03459
03460
03461 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if fill_i;
03462 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03463 {
03464
03465 LField<T,D> &fill = *(*fill_i).second;
03466 const NDIndex<D> &fill_alloc = fill.getAllocated();
03467 if ( slab.touches( fill_alloc ) )
03468 {
03469
03470 NDIndex<D> dest = slab.intersect( fill_alloc );
03471
03472
03473 NDIndex<D> src = dest;
03474 src[d] = offset - src[d];
03475
03476
03477 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if from_i;
03478 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03479 {
03480
03481 LField<T,D> &from = *(*from_i).second;
03482 const NDIndex<D> &from_owned = from.getOwned();
03483 const NDIndex<D> &from_alloc = from.getAllocated();
03484
03485 if ( src.touches( from_owned ) )
03486 {
03487
03488
03489
03490
03491
03492
03493
03494
03495
03496
03497
03498
03499 fill.Uncompress();
03500
03501 NDIndex<D> from_it = src.intersect( from_alloc );
03502 NDIndex<D> fill_it = dest.plugBase( from_it );
03503
03504 typedef typename LField<T,D>::iterator LFI;
03505 LFI lhs = fill.begin(fill_it);
03506 LFI rhs = from.begin(from_it);
03507
03508 BrickExpression<D,LFI,LFI,OpBCFunctionEq<T> >
03509 (lhs,rhs,OpBCFunctionEq<T>(ff.Func)).apply();
03510 }
03511 }
03512 }
03513 }
03514 }
03515
03517
03518
03519
03520
03521
03522
03523
03524 template<class T>
03525 struct OpBCFunctionEqComponent
03526 {
03527 OpBCFunctionEqComponent(typename ApplyToComponentType<T>::type (*func)
03528 ( typename ApplyToComponentType<T>::type ),
03529 int c) :
03530 Func(func), Component(c) {}
03531 typename ApplyToComponentType<T>::type
03532 (*Func)( typename ApplyToComponentType<T>::type );
03533 int Component;
03534 };
03535 template<class T>
03536 inline void PETE_apply(const OpBCFunctionEqComponent<T>& e, T& a, const T& b)
03537 { a[e.Component] = e.Func(b[e.Component]); }
03538
03539
03540
03541
03542
03543
03544
03545
03546 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,char)
03547 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,bool)
03548 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,int)
03549 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,unsigned)
03550 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,short)
03551 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,long)
03552 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,float)
03553 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,double)
03554 COMPONENT_APPLY_BUILTIN(OpBCFunctionEqComponent,dcomplex)
03555
03556
03558
03559
03560
03561
03562
03563
03564
03565
03566
03567
03568
03570
03571 template<class T, unsigned D, class M>
03572 void ComponentFunctionFaceBCApply(ComponentFunctionFace<T,D,M,Cell>& ff,
03573 Field<T,D,M,Cell>& A );
03574 template<class T, unsigned D, class M>
03575 void ComponentFunctionFaceBCApply(ComponentFunctionFace<T,D,M,Vert>& ff,
03576 Field<T,D,M,Vert>& A );
03577 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
03578 void ComponentFunctionFaceBCApply(ComponentFunctionFace<T,D,M,
03579 CartesianCentering<CE,D,NC> >& ff,
03580 Field<T,D,M,
03581 CartesianCentering<CE,D,NC> >& A );
03582
03583 template<class T, unsigned D, class M, class C>
03584 void ComponentFunctionFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
03585 {
03586 ComponentFunctionFaceBCApply(*this, A);
03587 }
03588
03589
03590
03591
03592
03593
03594 template<class T, unsigned D, class M>
03595 void ComponentFunctionFaceBCApply(ComponentFunctionFace<T,D,M,Cell>& ff,
03596 Field<T,D,M,Cell>& A )
03597 {
03598 TAU_TYPE_STRING(taustr, "void (" + CT(ff) + ", " + CT(A) + " )" );
03599 TAU_PROFILE("ComponentFunctionFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
03600
03601
03602
03603
03604
03605 const NDIndex<D>& domain( A.getDomain() );
03606 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
03607 unsigned d = ff.getFace()/2;
03608 int offset;
03609 if ( ff.getFace() & 1 ) {
03610
03611 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
03612 offset = 2*domain[d].max() + 1;
03613 }
03614 else {
03615 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03616 offset = 2*domain[d].min() - 1;
03617 }
03618
03619
03620 typename Field<T,D,M,Cell>::iterator_if fill_i;
03621 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03622 {
03623
03624 LField<T,D> &fill = *(*fill_i).second;
03625 const NDIndex<D> &fill_alloc = fill.getAllocated();
03626 if ( slab.touches( fill_alloc ) )
03627 {
03628
03629 NDIndex<D> dest = slab.intersect( fill_alloc );
03630
03631
03632 NDIndex<D> src = dest;
03633 src[d] = offset - src[d];
03634
03635
03636 typename Field<T,D,M,Cell>::iterator_if from_i;
03637 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03638 {
03639
03640 LField<T,D> &from = *(*from_i).second;
03641 const NDIndex<D> &from_owned = from.getOwned();
03642 const NDIndex<D> &from_alloc = from.getAllocated();
03643
03644 if ( src.touches( from_owned ) )
03645 {
03646
03647
03648
03649
03650
03651
03652
03653
03654
03655
03656
03657
03658 fill.Uncompress();
03659
03660 NDIndex<D> from_it = src.intersect( from_alloc );
03661 NDIndex<D> fill_it = dest.plugBase( from_it );
03662
03663 typedef typename LField<T,D>::iterator LFI;
03664 LFI lhs = fill.begin(fill_it);
03665 LFI rhs = from.begin(from_it);
03666
03667 BrickExpression<D,LFI,LFI,OpBCFunctionEqComponent<T> >
03668 (lhs,rhs,OpBCFunctionEqComponent<T>
03669 (ff.Func,ff.getComponent())).apply();
03670 }
03671 }
03672 }
03673 }
03674 }
03675
03676
03677
03678
03679
03680
03681 template<class T, unsigned D, class M>
03682 void ComponentFunctionFaceBCApply(ComponentFunctionFace<T,D,M,Vert>& ff,
03683 Field<T,D,M,Vert>& A )
03684 {
03685 TAU_TYPE_STRING(taustr, "void (" + CT(ff) + ", " + CT(A) + " )" );
03686 TAU_PROFILE("ComponentFunctionFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
03687
03688
03689
03690
03691
03692 const NDIndex<D>& domain( A.getDomain() );
03693 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
03694 unsigned d = ff.getFace()/2;
03695 int offset;
03696 if ( ff.getFace() & 1 ) {
03697
03698 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
03699
03700
03701 offset = 2*domain[d].max() + 1 - 1;
03702 }
03703 else {
03704 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03705
03706
03707 offset = 2*domain[d].min() - 1 + 1;
03708 }
03709
03710
03711 typename Field<T,D,M,Vert>::iterator_if fill_i;
03712 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03713 {
03714
03715 LField<T,D> &fill = *(*fill_i).second;
03716 const NDIndex<D> &fill_alloc = fill.getAllocated();
03717 if ( slab.touches( fill_alloc ) )
03718 {
03719
03720 NDIndex<D> dest = slab.intersect( fill_alloc );
03721
03722
03723 NDIndex<D> src = dest;
03724 src[d] = offset - src[d];
03725
03726
03727 typename Field<T,D,M,Vert>::iterator_if from_i;
03728 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03729 {
03730
03731 LField<T,D> &from = *(*from_i).second;
03732 const NDIndex<D> &from_owned = from.getOwned();
03733 const NDIndex<D> &from_alloc = from.getAllocated();
03734
03735 if ( src.touches( from_owned ) )
03736 {
03737
03738
03739
03740
03741
03742
03743
03744
03745
03746
03747
03748
03749 fill.Uncompress();
03750
03751 NDIndex<D> from_it = src.intersect( from_alloc );
03752 NDIndex<D> fill_it = dest.plugBase( from_it );
03753
03754 typedef typename LField<T,D>::iterator LFI;
03755 LFI lhs = fill.begin(fill_it);
03756 LFI rhs = from.begin(from_it);
03757
03758 BrickExpression<D,LFI,LFI,OpBCFunctionEqComponent<T> >
03759 (lhs,rhs,OpBCFunctionEqComponent<T>
03760 (ff.Func,ff.getComponent())).apply();
03761 }
03762 }
03763 }
03764 }
03765 }
03766
03767
03768
03769
03770
03771
03772 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
03773 void ComponentFunctionFaceBCApply(ComponentFunctionFace<T,D,M,
03774 CartesianCentering<CE,D,NC> >& ff,
03775 Field<T,D,M,CartesianCentering<CE,D,NC> >& A )
03776 {
03777 TAU_TYPE_STRING(taustr, "void (" + CT(ff) + ", " + CT(A) + " )" );
03778 TAU_PROFILE("ComponentFunctionFaceBCApply()", taustr, TAU_FIELD | TAU_ASSIGN);
03779
03780
03781
03782
03783
03784 const NDIndex<D>& domain( A.getDomain() );
03785 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
03786 unsigned d = ff.getFace()/2;
03787 int offset;
03788 if ( ff.getFace() & 1 ) {
03789
03790 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
03791
03792
03793 if (CE[ff.getComponent() + d*NC] == CELL) {
03794 ERRORMSG("check src code, had to change ef (not known) to ff ??? " << endl);
03795 offset = 2*domain[d].max() + 1;
03796 } else {
03797 offset = 2*domain[d].max() + 1 - 1;
03798 }
03799 }
03800 else {
03801 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
03802
03803
03804 if (CE[ff.getComponent() + d*NC] == CELL) {
03805 offset = 2*domain[d].min() - 1;
03806 } else {
03807 offset = 2*domain[d].min() - 1 + 1;
03808 }
03809 }
03810
03811
03812 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if fill_i;
03813 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i)
03814 {
03815
03816 LField<T,D> &fill = *(*fill_i).second;
03817 const NDIndex<D> &fill_alloc = fill.getAllocated();
03818 if ( slab.touches( fill_alloc ) )
03819 {
03820
03821 NDIndex<D> dest = slab.intersect( fill_alloc );
03822
03823
03824 NDIndex<D> src = dest;
03825 src[d] = offset - src[d];
03826
03827
03828 typename Field<T,D,M,CartesianCentering<CE,D,NC> >::iterator_if from_i;
03829 for (from_i=A.begin_if(); from_i!=A.end_if(); ++from_i)
03830 {
03831
03832 LField<T,D> &from = *(*from_i).second;
03833 const NDIndex<D> &from_owned = from.getOwned();
03834 const NDIndex<D> &from_alloc = from.getAllocated();
03835
03836 if ( src.touches( from_owned ) )
03837 {
03838
03839
03840
03841
03842
03843
03844
03845
03846
03847
03848
03849
03850 fill.Uncompress();
03851
03852 NDIndex<D> from_it = src.intersect( from_alloc );
03853 NDIndex<D> fill_it = dest.plugBase( from_it );
03854
03855 typedef typename LField<T,D>::iterator LFI;
03856 LFI lhs = fill.begin(fill_it);
03857 LFI rhs = from.begin(from_it);
03858
03859 BrickExpression<D,LFI,LFI,OpBCFunctionEqComponent<T> >
03860 (lhs,rhs,OpBCFunctionEqComponent<T>
03861 (ff.Func,ff.getComponent())).apply();
03862 }
03863 }
03864 }
03865 }
03866 }
03867
03869
03870
03871
03872
03873
03874
03875
03876
03877
03878
03879
03880
03881
03882
03883
03884
03885
03886 template<class T, unsigned int D, class M, class C>
03887 void EurekaFace<T,D,M,C>::apply(Field<T,D,M,C>& field)
03888 {
03889
03890
03891 NDIndex<D> slab = calcEurekaSlabToFill(field,BCondBase<T,D,M,C>::m_face,this->getComponent());
03892
03893
03894 fillSlabWithZero(field,slab,this->getComponent());
03895 }
03896
03898
03899
03900
03901
03902
03903
03904
03905
03906
03907
03908
03909
03910
03911
03912
03913
03914
03915 #ifdef __MWERKS__
03916
03917
03918
03919 template<class T, int D>
03920 #else
03921 template<class T>
03922 #endif // __MWERKS__
03923 struct EurekaAssign
03924 {
03925 static const T& get(const T& x, int) {
03926 ERRORMSG("Eureka assign to a component of class without an op[]!"<<endl);
03927 return x;
03928 }
03929 static T& get(T& x, int) {
03930 ERRORMSG("Eureka assign to a component of class without an op[]!"<<endl);
03931 return x;
03932 }
03933 int m_component;
03934 EurekaAssign(int c) : m_component(c) {}
03935 };
03936
03937
03938
03939
03940
03941
03942
03943
03944
03945 template<class T, unsigned int D, class M, class C>
03946 static void
03947 fillSlabWithZero(Field<T,D,M,C>& field,
03948 const NDIndex<D>& fillDomain,
03949 int component)
03950 {
03951
03952
03953
03954 typename Field<T,D,M,C>::iterator_if lp = field.begin_if();
03955 typename Field<T,D,M,C>::iterator_if done = field.end_if();
03956 for (; lp != done ; ++lp )
03957 {
03958
03959 LField<T,D>& lf = *(*lp).second;
03960
03961
03962 NDIndex<D> localDomain = lf.getAllocated();
03963
03964
03965 if ( fillDomain.touches(localDomain) )
03966 {
03967
03968 if ( lf.IsCompressed() )
03969 {
03970
03971 if ( component == BCondBase<T,D,M,C>::allComponents )
03972 {
03973
03974 if ( *lf.begin() == 0 )
03975 continue;
03976 }
03977
03978 else
03979 {
03980
03981
03982
03983 #ifdef __MWERKS__
03984
03985
03986
03987 if ( EurekaAssign<T,D>::get(*lf.begin(),component)==0 )
03988 continue;
03989 #else
03990 if ( EurekaAssign<T>::get(*lf.begin(),component)==0 )
03991 continue;
03992 #endif // __MWERKS__
03993 }
03994
03995
03996 lf.Uncompress();
03997 }
03998
03999
04000 NDIndex<D> intersectDomain = fillDomain.intersect(localDomain);
04001
04002
04003 typename LField<T,D>::iterator data = lf.begin(intersectDomain);
04004
04005
04006
04007
04008 typedef typename LField<T,D>::iterator Lhs_t;
04009 typedef PETE_Scalar<T> Rhs_t;
04010
04011
04012 if ( component == BCondBase<T,D,M,C>::allComponents )
04013 {
04014
04015 typedef BrickExpression<D,Lhs_t,Rhs_t,OpAssign> Expr_t;
04016
04017
04018 Expr_t(data,Rhs_t(0)).apply();
04019 }
04020
04021 else
04022 {
04023
04024 #ifdef __MWERKS__
04025
04026
04027
04028 typedef BrickExpression<D,Lhs_t,Rhs_t,EurekaAssign<T,D> > Expr_t;
04029 #else
04030 typedef BrickExpression<D,Lhs_t,Rhs_t,EurekaAssign<T> > Expr_t;
04031 #endif // __MWERKS__
04032
04033
04034 PAssert(component>=0);
04035
04036
04037 Expr_t(data,Rhs_t(0),component).apply();
04038
04039
04040 }
04041 }
04042 }
04043 }
04044
04045
04046
04047
04048
04049 template<class T, unsigned int D>
04050 #ifdef __MWERKS__
04051
04052
04053
04054 struct EurekaAssign< Vektor<T,D>, D >
04055 #else
04056 struct EurekaAssign< Vektor<T,D> >
04057 #endif
04058 {
04059 static T& get( Vektor<T,D>& x, int c ) { return x[c]; }
04060 static T get( const Vektor<T,D>& x, int c ) { return x[c]; }
04061 EurekaAssign(int c) : m_component(c) {}
04062 int m_component;
04063 };
04064
04065 template<class T, unsigned int D>
04066 #ifdef __MWERKS__
04067
04068
04069
04070 struct EurekaAssign< Tenzor<T,D>, D >
04071 #else
04072 struct EurekaAssign< Tenzor<T,D> >
04073 #endif
04074 {
04075 static T& get( Tenzor<T,D>& x, int c ) { return x[c]; }
04076 static T get( const Tenzor<T,D>& x, int c ) { return x[c]; }
04077 EurekaAssign(int c) : m_component(c) {}
04078 int m_component;
04079 };
04080
04081 template<class T, unsigned int D>
04082 #ifdef __MWERKS__
04083
04084
04085
04086 struct EurekaAssign< AntiSymTenzor<T,D>, D >
04087 #else
04088 struct EurekaAssign< AntiSymTenzor<T,D> >
04089 #endif
04090 {
04091 static T& get( AntiSymTenzor<T,D>& x, int c ) { return x[c]; }
04092 static T get( const AntiSymTenzor<T,D>& x, int c ) { return x[c]; }
04093 EurekaAssign(int c) : m_component(c) {}
04094 int m_component;
04095 };
04096
04097 template<class T, unsigned int D>
04098 #ifdef __MWERKS__
04099
04100
04101
04102 struct EurekaAssign< SymTenzor<T,D>, D >
04103 #else
04104 struct EurekaAssign< SymTenzor<T,D> >
04105 #endif
04106 {
04107 static T& get( SymTenzor<T,D>& x, int c ) { return x[c]; }
04108 static T get( const SymTenzor<T,D>& x, int c ) { return x[c]; }
04109 EurekaAssign(int c) : m_component(c) {}
04110 int m_component;
04111 };
04112
04113
04114
04115
04116
04117
04118 #ifdef __MWERKS__
04119
04120
04121
04122 template<class T, int D>
04123 inline void PETE_apply(const EurekaAssign<T,D>& e, T& a, const T& b)
04124 {
04125 EurekaAssign<T,D>::get(a,e.m_component) =
04126 EurekaAssign<T,D>::get(b,e.m_component);
04127 }
04128 #else
04129 template<class T>
04130 inline void PETE_apply(const EurekaAssign<T>& e, T& a, const T& b)
04131 {
04132 EurekaAssign<T>::get(a,e.m_component) =
04133 EurekaAssign<T>::get(b,e.m_component);
04134 }
04135 #endif // __MWERKS__
04136
04137
04139
04140
04141
04142
04143
04144
04145
04146 template<unsigned int D>
04147 inline NDIndex<D>
04148 calcEurekaDomain(const NDIndex<D>& realDomain,
04149 int face,
04150 const GuardCellSizes<D>& gc)
04151 {
04152 NDIndex<D> slab = AddGuardCells(realDomain,gc);
04153
04154
04155 int dim = face/2;
04156
04157
04158 int low,high;
04159
04160
04161 if ( face&1 )
04162 {
04163
04164 high = slab[dim].max();
04165 low = high - gc.right(dim);
04166 }
04167
04168 else
04169 {
04170
04171 low = slab[dim].min();
04172 high = low + gc.left(dim);
04173 }
04174
04175
04176 PAssert( low<=high );
04177 PAssert( high<=slab[dim].max() );
04178 PAssert( low >=slab[dim].min() );
04179
04180
04181 slab[dim] = Index(low,high);
04182 return slab;
04183 }
04184
04186
04187
04188
04189
04190
04191
04192
04193 template<class T, unsigned int D, class M>
04194 static NDIndex<D>
04195 calcEurekaSlabToFill(const Field<T,D,M,Cell>& field, int face,int)
04196 {
04197 return calcEurekaDomain(field.getDomain(),face,field.getGC());
04198 }
04199
04200 template<class T, unsigned int D, class M>
04201 static NDIndex<D>
04202 calcEurekaSlabToFill(const Field<T,D,M,Vert>& field, int face,int)
04203 {
04204 return calcEurekaDomain(field.getDomain(),face,field.getGC());
04205 }
04206
04208
04209
04210
04211
04212
04213
04214
04215
04216
04217
04218
04219 template<class T, unsigned D, class M, const CenteringEnum* CE, unsigned NC>
04220 static NDIndex<D>
04221 calcEurekaSlabToFill(const Field<T,D,M,CartesianCentering<CE,D,NC> >& field,
04222 int face,
04223 int component)
04224 {
04225
04226 typedef CartesianCentering<CE,D,NC> C;
04227
04228
04229 int d = face/2;
04230
04231
04232 NDIndex<D> slab;
04233
04234
04235 if ( component == BCondBase<T,D,M,C>::allComponents )
04236 {
04237
04238 slab = calcEurekaDomain(field.getDomain(),face,field.getGC());
04239 }
04240
04241 else
04242 {
04243
04244
04245 slab = AddGuardCells(field.getDomain(),field.getGC());
04246
04247
04248 if ( face&1 )
04249 {
04250
04251 int low = field.getDomain()[d].min() + field.get_mesh().gridSizes[d] - 1;
04252 int high = low + field.rightGuard(d);
04253
04254
04255 if (CE[component + d*NC] == CELL)
04256 low++;
04257
04258
04259 PAssert( low<=high );
04260 PAssert( high<=slab[d].max() );
04261 PAssert( low >=slab[d].min() );
04262
04263
04264 slab[d] = Index(low,high);
04265 }
04266
04267 else
04268 {
04269
04270 int low = slab[d].min();
04271 int high = low + field.leftGuard(d) ;
04272
04273
04274 if (CE[component + d*NC] == CELL)
04275 high--;
04276
04277
04278 PAssert( low<=high );
04279 PAssert( high<=slab[d].max() );
04280 PAssert( low >=slab[d].min() );
04281
04282
04283 slab[d] = Index(low,high);
04284 }
04285 }
04286
04287
04288 return slab;
04289 }
04290
04292
04293
04294
04295
04296
04297
04298
04299
04300
04301
04302
04303
04304
04305
04306
04308
04309 template<class T, unsigned D, class M, class C>
04310 void LinearExtrapolateFaceBCApply(LinearExtrapolateFace<T,D,M,C>& ef,
04311 Field<T,D,M,C>& A );
04312
04313 template<class T, unsigned D, class M, class C>
04314 void LinearExtrapolateFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
04315 {
04316 LinearExtrapolateFaceBCApply(*this, A);
04317 }
04318
04319
04320 template<class T, unsigned D, class M, class C>
04321 void LinearExtrapolateFaceBCApply2(const NDIndex<D> &dest,
04322 const NDIndex<D> &src1,
04323 const NDIndex<D> &src2,
04324 LField<T,D> &fill,
04325 LinearExtrapolateFace<T,D,M,C> &ef,
04326 int slopeMultipplier)
04327 {
04328
04329
04330
04331
04332
04333 if (fill.IsCompressed()) { return; }
04334
04335
04336 typedef typename LField<T,D>::iterator LFI;
04337 LFI lhs = fill.begin(dest);
04338 LFI rhs1 = fill.begin(src1);
04339 LFI rhs2 = fill.begin(src2);
04340 LFI endi = fill.end();
04341
04342
04343
04344
04345 for ( ; lhs != endi, rhs1 != endi, rhs2 != endi;
04346 ++lhs, ++rhs1, ++rhs2) {
04347 *lhs = (*rhs2 - *rhs1)*slopeMultipplier + *rhs1;
04348 }
04349
04350 }
04351
04352
04353
04354
04355
04356
04357
04358
04359 template<class T, unsigned D, class M, class C>
04360 void LinearExtrapolateFaceBCApply(LinearExtrapolateFace<T,D,M,C>& ef,
04361 Field<T,D,M,C>& A )
04362 {
04363 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
04364 TAU_PROFILE("LinearExtrapolateFaceBCApply()", taustr,TAU_FIELD | TAU_ASSIGN);
04365
04366
04367
04368
04369
04370 const NDIndex<D>& domain( A.getDomain() );
04371 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
04372
04373
04374
04375
04376
04377 unsigned d = ef.getFace()/2;
04378
04379
04380 int nGuardLayers;
04381
04382
04383
04384
04385
04386
04387 if (ef.getFace() & 1) {
04388
04389
04390
04391 nGuardLayers = A.rightGuard(d);
04392
04393 } else {
04394
04395
04396
04397
04398 nGuardLayers = A.leftGuard(d);
04399
04400 }
04401
04402
04403
04404
04405 for (int guardLayer = 1; guardLayer <= nGuardLayers; guardLayer++) {
04406
04407
04408
04409 int slopeMultipplier = -1*guardLayer;
04410
04411 if (ef.getFace() & 1) {
04412 slab[d] = Index(domain[d].max() + guardLayer,
04413 domain[d].max() + guardLayer);
04414 } else {
04415 slab[d] = Index(domain[d].min() - guardLayer,
04416 domain[d].min() - guardLayer);
04417 }
04418
04419
04420
04421 typename Field<T,D,M,Cell>::iterator_if fill_i;
04422 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i) {
04423
04424
04425
04426
04427 LField<T,D> &fill = *(*fill_i).second;
04428
04429
04430 const NDIndex<D> &fill_alloc = fill.getAllocated();
04431
04432
04433
04434
04435
04436
04437 if (slab.touches(fill_alloc)) {
04438
04439
04440 NDIndex<D> dest = slab.intersect(fill_alloc);
04441
04442
04443
04444
04445
04446
04447
04448
04449 NDIndex<D> src1 = dest;
04450 NDIndex<D> src2 = dest;
04451 if (ef.getFace() & 1) {
04452 src2[d] = Index(domain[d].max() - 1, domain[d].max() - 1, 1);
04453 src1[d] = Index(domain[d].max(), domain[d].max(), 1);
04454 } else {
04455 src1[d] = Index(0,0,1);
04456 src2[d] = Index(1,1,1);
04457 }
04458
04459
04460
04461
04462
04463 LinearExtrapolateFaceBCApply2(dest, src1, src2, fill, ef,
04464 slopeMultipplier);
04465
04466 }
04467 }
04468 }
04469 }
04470
04471
04473
04474
04475
04476
04477
04478
04479
04480
04481
04482
04483
04484
04485
04486
04487
04488
04489 template<class T, unsigned D, class M, class C>
04490 void ComponentLinearExtrapolateFaceBCApply(ComponentLinearExtrapolateFace
04491 <T,D,M,C>& ef,
04492 Field<T,D,M,C>& A );
04493
04494 template<class T, unsigned D, class M, class C>
04495 void ComponentLinearExtrapolateFace<T,D,M,C>::apply( Field<T,D,M,C>& A )
04496 {
04497 ComponentLinearExtrapolateFaceBCApply(*this, A);
04498 }
04499
04500
04501 template<class T, unsigned D, class M, class C>
04502 void ComponentLinearExtrapolateFaceBCApply2(const NDIndex<D> &dest,
04503 const NDIndex<D> &src1,
04504 const NDIndex<D> &src2,
04505 LField<T,D> &fill,
04506 ComponentLinearExtrapolateFace
04507 <T,D,M,C> &ef,
04508 int slopeMultipplier)
04509 {
04510
04511
04512
04513
04514
04515 if (fill.IsCompressed()) { return; }
04516
04517
04518 typedef typename LField<T,D>::iterator LFI;
04519 LFI lhs = fill.begin(dest);
04520 LFI rhs1 = fill.begin(src1);
04521 LFI rhs2 = fill.begin(src2);
04522 LFI endi = fill.end();
04523
04524
04525
04526
04527
04528 int component = ef.getComponent();
04529 for ( ; lhs != endi, rhs1 != endi, rhs2 != endi;
04530 ++lhs, ++rhs1, ++rhs2) {
04531 (*lhs)[component] =
04532 ((*rhs2)[component] - (*rhs1)[component])*slopeMultipplier +
04533 (*rhs1)[component];
04534 }
04535
04536 }
04537
04538
04539
04540
04541
04542
04543
04544
04545 template<class T, unsigned D, class M, class C>
04546 void ComponentLinearExtrapolateFaceBCApply(ComponentLinearExtrapolateFace
04547 <T,D,M,C>& ef,
04548 Field<T,D,M,C>& A )
04549 {
04550 TAU_TYPE_STRING(taustr, "void (" + CT(ef) + ", " + CT(A) + " )" );
04551 TAU_PROFILE("ComponentLinearExtrapolateFaceBCApply()",
04552 taustr,TAU_FIELD | TAU_ASSIGN);
04553
04554
04555
04556
04557
04558 const NDIndex<D>& domain( A.getDomain() );
04559 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
04560
04561
04562
04563
04564
04565 unsigned d = ef.getFace()/2;
04566
04567
04568 int nGuardLayers;
04569
04570
04571
04572
04573
04574
04575 if (ef.getFace() & 1) {
04576
04577
04578
04579 nGuardLayers = A.rightGuard(d);
04580
04581 } else {
04582
04583
04584
04585
04586 nGuardLayers = A.leftGuard(d);
04587
04588 }
04589
04590
04591
04592
04593 for (int guardLayer = 1; guardLayer <= nGuardLayers; guardLayer++) {
04594
04595
04596
04597 int slopeMultipplier = -1*guardLayer;
04598
04599 if (ef.getFace() & 1) {
04600 slab[d] = Index(domain[d].max() + guardLayer,
04601 domain[d].max() + guardLayer);
04602 } else {
04603 slab[d] = Index(domain[d].min() - guardLayer,
04604 domain[d].min() - guardLayer);
04605 }
04606
04607
04608
04609 typename Field<T,D,M,Cell>::iterator_if fill_i;
04610 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i) {
04611
04612
04613
04614
04615 LField<T,D> &fill = *(*fill_i).second;
04616
04617
04618 const NDIndex<D> &fill_alloc = fill.getAllocated();
04619
04620
04621
04622
04623
04624
04625 if (slab.touches(fill_alloc)) {
04626
04627
04628 NDIndex<D> dest = slab.intersect(fill_alloc);
04629
04630
04631
04632
04633
04634
04635
04636
04637 NDIndex<D> src1 = dest;
04638 NDIndex<D> src2 = dest;
04639 if (ef.getFace() & 1) {
04640 src2[d] = Index(domain[d].max() - 1, domain[d].max() - 1, 1);
04641 src1[d] = Index(domain[d].max(), domain[d].max(), 1);
04642 } else {
04643 src1[d] = Index(0,0,1);
04644 src2[d] = Index(1,1,1);
04645 }
04646
04647
04648
04649
04650
04651 ComponentLinearExtrapolateFaceBCApply2(dest, src1, src2, fill, ef,
04652 slopeMultipplier);
04653
04654 }
04655 }
04656 }
04657 }
04658
04659
04660
04661
04662 template<class T, unsigned D, class M, class C>
04663 PatchBC<T,D,M,C>::
04664 PatchBC(unsigned face)
04665 : BCondBase<T,D,M,C>(face)
04666 {
04667 TAU_TYPE_STRING(taustr, "void ( unsigned )" );
04668 TAU_PROFILE("PatchBC::PatchBC()", taustr, TAU_FIELD);
04669 }
04670
04671
04672
04673
04674
04675
04676
04677
04678 template<class T, unsigned D, class M, class C>
04679 void PatchBC<T,D,M,C>::apply( Field<T,D,M,C>& A )
04680 {
04681 TAU_TYPE_STRING(taustr, "void (" + CT(A) + " )" );
04682 TAU_PROFILE("PatchBC()", taustr, TAU_FIELD | TAU_ASSIGN);
04683
04684
04685
04686
04687
04688
04689
04690
04691 const NDIndex<D>& domain( A.getDomain() );
04692
04693
04694
04695
04696
04697 NDIndex<D> slab = AddGuardCells(domain,A.getGuardCellSizes());
04698
04699
04700
04701
04702
04703 unsigned d = this->getFace()/2;
04704
04705
04706
04707
04708
04709
04710 if ( this->getFace() & 1 )
04711 slab[d] = Index( domain[d].max() + 1, domain[d].max() + A.rightGuard(d));
04712 else
04713 slab[d] = Index( domain[d].min() - A.leftGuard(d), domain[d].min()-1 );
04714
04715
04716
04717
04718 int lfindex = 0;
04719 typename Field<T,D,M,C>::iterator_if fill_i;
04720 typename Field<T,D,M,C>::iterator p = A.begin();
04721 for (fill_i=A.begin_if(); fill_i!=A.end_if(); ++fill_i, ++lfindex)
04722 {
04723
04724
04725
04726
04727
04728 LField<T,D> &fill = *(*fill_i).second;
04729 const NDIndex<D> &fill_alloc = fill.getAllocated();
04730 const NDIndex<D> &fill_owned = fill.getOwned();
04731
04732
04733
04734
04735 if ( slab.touches( fill_alloc ) )
04736 {
04737
04738 fill.Uncompress();
04739
04740
04741
04742
04743 NDIndex<D> dest = slab.intersect( fill_alloc );
04744
04745
04746
04747
04748 vec<int,D> v;
04749 for (int i=0; i<D; ++i)
04750 v[i] = dest[i].first();
04751 p.SetCurrentLocation( FieldLoc<D>(v,lfindex) );
04752
04753
04754
04755
04756 applyPatch(p,dest);
04757 }
04758 }
04759 }
04760
04761
04762
04763 #undef COMPONENT_APPLY_BUILTIN
04764
04765
04766
04767
04768
04769