-
Notifications
You must be signed in to change notification settings - Fork 0
/
MaterialRunlineEncoder.cpp
154 lines (137 loc) · 5.36 KB
/
MaterialRunlineEncoder.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
/*
* MaterialSetupUpdateEquation.cpp
* TROGDOR
*
* Created by Paul Hansen on 8/19/09.
* Copyright 2009 Stanford University. All rights reserved.
*
*/
#include "MaterialRunlineEncoder.h"
#include "VoxelizedPartition.h"
#include "YeeUtilities.h"
#include <iostream>
using namespace std;
using namespace YeeUtilities;
BulkMaterialRLE::
BulkMaterialRLE()
{
setNeighborFieldContinuity(kNeighborFieldTransverse4);
}
void BulkMaterialRLE::
endRunline(const VoxelizedPartition & vp, const Vector3i & lastHalfCell)
{
// LOG << "Ending runline from " << firstHalfCell() << " length " << length()
// << endl;
Paint* paint = vp.voxels()(firstHalfCell());
int oct = octant(firstHalfCell());
int dir0 = xyz(oct); // this is the field direction, "i"
int dir1 = (dir0+1)%3; // first transverse direction, "j"
int dir2 = (dir1+1)%3; // second transverse direction, "k"
SBMRunline* rl = new SBMRunline;
rl->length = length();
rl->auxIndex = vp.indices()(firstHalfCell());
rl->f_i = vp.lattice().pointer(firstHalfCell());
// WATCH CAREFULLY
//
// f_i will be, e.g., Ex
// f_j will be Hy
// f_k will be Hz
//
// The non-trivial piece here is that the neighbor in direction j is a
// field that points along k. I've screwed this up several times!
if (paint->hasCurlBuffer(2*dir2))
rl->f_j[0] = paint->curlBuffer(2*dir2)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir2));
else
rl->f_j[0] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir2));
if (paint->hasCurlBuffer(2*dir2+1))
rl->f_j[1] = paint->curlBuffer(2*dir2+1)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir2+1));
else
rl->f_j[1] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir2+1));
if (paint->hasCurlBuffer(2*dir1))
rl->f_k[0] = paint->curlBuffer(2*dir1)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir1));
else
rl->f_k[0] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir1));
if (paint->hasCurlBuffer(2*dir1+1))
rl->f_k[1] = paint->curlBuffer(2*dir1+1)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir1+1));
else
rl->f_k[1] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir1+1));
if (isE(oct))
mRunlinesE[dir0].push_back(SBMRunlinePtr(rl));
else
mRunlinesH[dir0].push_back(SBMRunlinePtr(rl));
}
BulkPMLRLE::
BulkPMLRLE()
{
setNeighborFieldContinuity(kNeighborFieldTransverse4);
setLineContinuity(kLineContinuityRequired);
}
void BulkPMLRLE::
endRunline(const VoxelizedPartition & vp, const Vector3i & lastHalfCell)
{
Paint* paint = vp.voxels()(firstHalfCell());
int oct = octant(firstHalfCell());
int dir0 = xyz(oct); // this is the field direction, "i"
int dir1 = (dir0+1)%3; // first transverse direction, "j"
int dir2 = (dir1+1)%3; // second transverse direction, "k"
SBPMRunline* rl = new SBPMRunline;
rl->length = length();
rl->auxIndex = vp.indices()(firstHalfCell());
rl->f_i = vp.lattice().pointer(firstHalfCell());
// WATCH CAREFULLY
//
// f_i will be, e.g., Ex
// f_j will be Hy
// f_k will be Hz
//
// The non-trivial piece here is that the neighbor in direction j is a
// field that points along k. I've screwed this up several times!
if (paint->hasCurlBuffer(2*dir2))
rl->f_j[0] = paint->curlBuffer(2*dir2)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir2));
else
rl->f_j[0] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir2));
if (paint->hasCurlBuffer(2*dir2+1))
rl->f_j[1] = paint->curlBuffer(2*dir2+1)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir2+1));
else
rl->f_j[1] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir2+1));
if (paint->hasCurlBuffer(2*dir1))
rl->f_k[0] = paint->curlBuffer(2*dir1)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir1));
else
rl->f_k[0] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir1));
if (paint->hasCurlBuffer(2*dir1+1))
rl->f_k[1] = paint->curlBuffer(2*dir1+1)->lattice()->wrappedPointer(
firstHalfCell()+cardinal(2*dir1+1));
else
rl->f_k[1] = vp.lattice().wrappedPointer(
firstHalfCell()+cardinal(2*dir1+1));
// PML aux stuff
// The start point of the runline *may* be outside the grid, *if* we are
// performing calculations on ghost points! This may happen in data-push
// adjoint update equations. In any case, usually the wrap does nothing.
Vector3i pmlDirections = vp.voxels()(firstHalfCell())->pmlDirections();
assert(vec_eq(vp.gridHalfCells().p1, 0));
Vector3i gridNumHalfCells = vp.gridHalfCells().size() + 1;
Vector3i wrappedStart = (firstHalfCell()+gridNumHalfCells)
% gridNumHalfCells;
assert(octant(wrappedStart) == octant(firstHalfCell()));
Rect3i pmlYeeCells = halfToYee(vp.pmlHalfCells(pmlDirections), oct);
rl->pmlDepthIndex = halfToYee(wrappedStart) - pmlYeeCells.p1;
if (isE(oct))
mRunlinesE[dir0].push_back(SBPMRunlinePtr(rl));
else
mRunlinesH[dir0].push_back(SBPMRunlinePtr(rl));
}