1 /*
2 *	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 *	This program is free software: you can redistribute it and/or modify
5 *	it under the terms of the GNU General Public License as published by
6 *	the Free Software Foundation, either version 3 of the License, or
7 *	(at your option) any later version.
8 *
9 *	This program is distributed in the hope that it will be useful,
10 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
11 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 *	GNU General Public License for more details.
13 *
14 *	You should have received a copy of the GNU General Public License
15 *	along with this program.  If not, see <http://www.gnu.org/licenses/>.
16 */
17 
18 #include "AdrOp.h"
19 
20 using namespace std;
21 
AdrOp(unsigned int muiImax,unsigned int muiJmax,unsigned int muiKmax,unsigned int muiLmax)22 AdrOp::AdrOp(unsigned int muiImax, unsigned int muiJmax, unsigned int muiKmax, unsigned int muiLmax)
23 {
24 	//error-handling...
25 	error = new ErrorMsg(9);
26 	if (error==NULL)
27 	{
28 		fprintf(stderr,"Memory allocation failed!! exiting...");
29 		exit(1);
30 	}
31 	error->SetMsg(1,"Adress Operator: Memory allocation failed!! exiting...");
32 	error->SetMsg(2,"Adress Operator: Invalid Adress requested!! exiting...");
33 	error->SetMsg(3,"Adress Operator: Invalid Position set!! exiting...");
34 	error->SetMsg(4,"Adress Operator: Invalid jump or passing end of iteration!! exiting...");
35 	error->SetMsg(5,"Adress Operator: 4D not yet implemented!! exiting...");
36 	error->SetMsg(6,"Adress Operator: Position not set!! exiting...");
37 	error->SetMsg(7,"Adress Operator: Cells not added to Adress Operator!! exiting...");
38 	error->SetMsg(8,"Adress Operator: Invalid Node!! exiting...");
39 	error->SetMsg(9,"Adress Operator: Grid invalid!! exiting...");
40 
41 	//if (muiImax<0) muiImax=0;
42 	//if (muiJmax<0) muiJmax=0;
43 	//if (muiKmax<0) muiKmax=0;
44 	//if (muiLmax<0) muiLmax=0;
45 
46 	uiDimension=0;
47 	if (muiImax>0) uiDimension++;
48 	else exit(-1);
49 	if (muiJmax>0) uiDimension++;
50 	else exit(-2);
51 	if (muiKmax>0) uiDimension++;
52 	if ( (muiLmax>0) && (muiKmax>0) ) uiDimension++;
53 //	cout << "\n-----Adress Operator created: Dimension: " << uiDimension << "----" <<endl;
54 	uiImax=muiImax;
55 	uiJmax=muiJmax;
56 	uiKmax=muiKmax;
57 	uiLmax=muiLmax=0;
58 	if (uiDimension==2) uiSize=uiImax*uiJmax;
59 	else if (uiDimension==3) uiSize=uiImax*uiJmax*uiKmax;
60 	else if (uiDimension==4) uiSize=uiImax*uiJmax*uiKmax*uiLmax;
61 	else uiSize=0;
62 	bPosSet=false;
63 	if (uiDimension==4) error->Error(5);
64 	iIshift=iJshift=iKshift=0;
65 	reflect=false;
66 	uiTypeOffset=0;
67 	clCellAdr=NULL;
68 	dGrid[0]=NULL;
69 	dGrid[1]=NULL;
70 	dGrid[2]=NULL;
71 	dGrid[3]=NULL;
72 	dDeltaUnit=1;
73 	bDebug=false;
74 }
75 
AdrOp(AdrOp * origOP)76 AdrOp::AdrOp(AdrOp* origOP)
77 {
78 	clCellAdr=NULL;
79 	error=NULL; // has to be done!!!
80 
81 	uiDimension=origOP->uiDimension;
82 	uiSize=origOP->uiSize;
83 	uiImax=origOP->uiImax;
84 	uiJmax=origOP->uiJmax;
85 	uiKmax=origOP->uiKmax;
86 	uiLmax=origOP->uiLmax;
87 	uiIpos=origOP->uiIpos;
88 	uiJpos=origOP->uiJpos;
89 	uiKpos=origOP->uiKpos;
90 	uiLpos=origOP->uiLpos;
91 	for (int ii=0; ii<4; ++ii) dGrid[ii]=origOP->dGrid[ii];
92 	dDeltaUnit=origOP->dDeltaUnit;
93 	iIshift=origOP->iIshift;
94 	iJshift=origOP->iJshift;
95 	iKshift=origOP->iKshift;
96 	for (int ii=0; ii<3; ++ii) iCellShift[ii]=origOP->iCellShift[ii];
97 	i=origOP->i;
98 	j=origOP->j;
99 	k=origOP->k;
100 	l=origOP->l;
101 	reflect=origOP->reflect;
102 	uiTypeOffset=origOP->uiTypeOffset;
103 
104 	bPosSet=origOP->bPosSet;
105 	bDebug=origOP->bDebug;
106 //	return;
107 	if (origOP->clCellAdr!=NULL) clCellAdr= new AdrOp(origOP->clCellAdr);
108 }
109 
~AdrOp()110 AdrOp::~AdrOp()
111 {
112 //	cerr << "\n------Adress Operator deconstructed-----\n" << endl;
113 	delete error;
114 	error=NULL;
115 	delete clCellAdr;
116 	clCellAdr=NULL;
117 }
118 
SetPos(unsigned int muiIpos,unsigned int muiJpos,unsigned int muiKpos,unsigned int muiLpos)119 unsigned int AdrOp::SetPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
120 {
121 	if (bDebug) fprintf(stderr,"AdrOp Debug:: SetPos(%d,%d,%d,%d) Max(%d,%d,%d,%d) \n",muiIpos,muiJpos,muiKpos,muiLpos,uiImax,uiJmax,uiKmax,uiLmax);
122 	bPosSet=false;
123 	if (muiIpos<uiImax) uiIpos=muiIpos;
124 	else error->Error(3);
125 	if (muiJpos<uiJmax) uiJpos=muiJpos;
126 	else error->Error(3);
127 	if ((muiKpos>=uiKmax) && (uiDimension>2)) error->Error(3);
128 	else if (uiDimension>2) uiKpos=muiKpos;
129 	if ((muiLpos>=uiLmax) && (uiDimension>3)) error->Error(3);
130 	else if (uiDimension>3) uiLpos=muiLpos;
131 	bPosSet=true;
132 //	cerr << "Position i:" << uiIpos << " j: " << uiJpos << " k: " << uiKpos << " l: " << 0 << " MAX: i:" << uiImax << " j: "  << uiJmax << " k: " << uiKmax << endl; //debug
133 	ADRESSEXPENSE(0,0,0,0,uiDimension+1,18)
134 	return GetPos();
135 }
136 
SetPosChecked(unsigned int muiIpos,unsigned int muiJpos,unsigned int muiKpos,unsigned int muiLpos)137 bool  AdrOp::SetPosChecked(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
138 {
139 	bPosSet=true;
140 	if (muiIpos<uiImax) uiIpos=muiIpos;
141 	else bPosSet=false;
142 	if (muiJpos<uiJmax) uiJpos=muiJpos;
143 	else bPosSet=false;
144 	if ((muiKpos>=uiKmax) && (uiDimension>2)) bPosSet=false;
145 	else if (uiDimension>2) uiKpos=muiKpos;
146 	if ((muiLpos>=uiLmax) && (uiDimension>3)) bPosSet=false;
147 	else if (uiDimension>3) uiLpos=muiLpos;
148 	ADRESSEXPENSE(0,0,0,0,uiDimension+1,18)
149 	return bPosSet;
150 }
151 
SetGrid(double * gridI,double * gridJ,double * gridK,double * gridL)152 void AdrOp::SetGrid(double *gridI,double *gridJ,double *gridK,double *gridL)
153 {
154 	dGrid[0]=gridI;
155 	dGrid[1]=gridJ;
156 	dGrid[2]=gridK;
157 	dGrid[3]=gridL;
158 	ADRESSEXPENSE(0,0,0,0,4,0)
159 }
160 
CheckPos(unsigned int muiIpos,unsigned int muiJpos,unsigned int muiKpos,unsigned int muiLpos)161 bool AdrOp::CheckPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
162 {
163 	bPosSet=true;
164 	if ((muiIpos>=uiImax)) bPosSet=false;
165 	if ((muiJpos>=uiJmax)) bPosSet=false;
166 	if ((muiKpos>=uiKmax) && (uiDimension>2)) bPosSet=false;
167 	if ((muiLpos>=uiLmax) && (uiDimension>3)) bPosSet=false;
168 	ADRESSEXPENSE(0,0,0,0,uiDimension+1,18)
169 	return bPosSet;
170 }
171 
CheckRelativePos(int muiIrel,int muiJrel,int muiKrel,int muiLrel)172 bool AdrOp::CheckRelativePos(int muiIrel,int muiJrel,int muiKrel, int muiLrel)
173 {
174 	bPosSet=true;
175 	if ((muiIrel+(int)uiIpos<0) || (muiIrel+(int)uiIpos>=(int)uiImax)) bPosSet=false;
176 	if ((muiJrel+(int)uiJpos<0) || (muiJrel+(int)uiJpos>=(int)uiJmax)) bPosSet=false;
177 	if (((muiKrel+(int)uiKpos<0) || (muiKrel+(int)uiKpos>=(int)uiKmax)) && (uiDimension>2)) bPosSet=false;
178 	if (((muiLrel+(int)uiLpos<0) || (muiLrel+(int)uiLpos>=(int)uiLmax)) && (uiDimension>3)) bPosSet=false;
179 	ADRESSEXPENSE(2*uiDimension,0,0,0,uiDimension+1,18)
180 	return bPosSet;
181 }
182 
GetPos(int muiIrel,int muiJrel,int muiKrel,int)183 unsigned int AdrOp::GetPos(int muiIrel, int muiJrel, int muiKrel, int /*muiLrel*/)
184 {
185 	if (bPosSet==false) error->Error(6);
186 	if (reflect)
187 	{
188 #if EXPENSE_LOG==1
189 		if (muiIrel+(int)uiIpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
190 			if (muiIrel+(int)uiIpos>(int)uiImax-1) ADRESSEXPENSE(4,1,0,0,1,0)
191 				if (muiJrel+(int)uiJpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
192 					if (muiJrel+(int)uiJpos>(int)uiJmax-1) ADRESSEXPENSE(4,1,0,0,1,0)
193 						if (muiKrel+(int)uiKpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
194 							if (muiKrel+(int)uiKpos>(int)uiKmax-1) ADRESSEXPENSE(4,1,0,0,1,0)
195 #endif
196 
197 								if (muiIrel+(int)uiIpos<0) muiIrel=-2*uiIpos-muiIrel-uiTypeOffset;
198 		if (muiIrel+(int)uiIpos>(int)uiImax-1) muiIrel=2*(uiImax-1-uiIpos)-muiIrel+uiTypeOffset;
199 		if (muiJrel+(int)uiJpos<0) muiJrel=-2*uiJpos-muiJrel-uiTypeOffset;
200 		if (muiJrel+(int)uiJpos>(int)uiJmax-1) muiJrel=2*(uiJmax-1-uiJpos)-muiJrel+uiTypeOffset;
201 		if (muiKrel+(int)uiKpos<0) muiKrel=-2*uiKpos-muiKrel-uiTypeOffset;
202 		if (muiKrel+(int)uiKpos>(int)uiKmax-1) muiKrel=2*(uiKmax-1-uiKpos)-muiKrel+uiTypeOffset;
203 	}
204 	if (uiDimension==2)
205 	{
206 		ADRESSEXPENSE(7,1,0,0,0,7)
207 		if ( (muiIrel+uiIpos<uiImax) && (muiJrel+uiJpos<uiJmax) )
208 			return (muiIrel+uiIpos)+(muiJrel+uiJpos)*uiImax;
209 		else error->Error(2);
210 		return 0;
211 	}
212 	else if (uiDimension==3)
213 	{
214 		ADRESSEXPENSE(9,3,0,0,0,11)
215 		if ( (muiIrel+uiIpos<uiImax) && (muiJrel+uiJpos<uiJmax) && (muiKrel+uiKpos<uiKmax) )
216 			return (muiIrel+uiIpos) + (muiJrel+uiJpos)*uiImax + (muiKrel+uiKpos)*uiJmax*uiImax;
217 		else error->Error(2);
218 		return 0;
219 	}
220 	else return 0;
221 }
222 
GetPosFromNode(int ny,unsigned int uiNode)223 unsigned int AdrOp::GetPosFromNode(int ny, unsigned int uiNode)
224 {
225 	while (ny<0) ny+=uiDimension;
226 	ny=ny%uiDimension;
227 	unsigned int help=uiNode,i=0,j=0,k=0,l=0;
228 	i=help%uiImax;
229 	help=help/uiImax;
230 	j=help%uiJmax;
231 	help=help/uiJmax;
232 	if (uiKmax>0)
233 	{
234 		k=help%uiKmax;
235 		help=help/uiKmax;
236 		l=help;
237 	}
238 	if (!CheckPos(i,j,k,l)) error->Error(8);
239 	ADRESSEXPENSE(0,7,0,0,13,4)
240 	switch (ny)
241 	{
242 	case 0:
243 		{
244 			return i;
245 			break;
246 		}
247 	case 1:
248 		{
249 			return j;
250 			break;
251 		}
252 	case 2:
253 		{
254 			return k;
255 			break;
256 		}
257 	case 3:
258 		{
259 			return l;
260 			break;
261 		}
262 	}
263 	return 0;
264 }
265 
GetNodeVolume(unsigned int uiNode)266 double AdrOp::GetNodeVolume(unsigned int uiNode)
267 {
268 	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
269 	double dVol=1;
270 	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
271 	unsigned int uiPos[4]={GetPosFromNode(0,uiNode),GetPosFromNode(1,uiNode),GetPosFromNode(2,uiNode),GetPosFromNode(3,uiNode)};
272 	for (unsigned int n=0; n<uiDimension; n++)
273 	{
274 		if ((uiPos[n]>0) && (uiPos[n]<uiMax[n]-1))
275 		{
276 			dVol*=0.5*dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]-1]);
277 			ADRESSEXPENSE(4,0,1,3,0,4)
278 		}
279 		else if ((uiPos[n]==0) && (uiPos[n]<uiMax[n]-1))
280 		{
281 			dVol*=dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]]);
282 			ADRESSEXPENSE(3,0,1,2,0,4)
283 		}
284 		else if ((uiPos[n]>0) && (uiPos[n]==uiMax[n]-1))
285 		{
286 			dVol*=dDeltaUnit*(dGrid[n][uiPos[n]]-dGrid[n][uiPos[n]-1]);
287 			ADRESSEXPENSE(3,0,1,2,0,4)
288 		}
289 	}
290 	return dVol;
291 }
292 
GetIndexWidth(int ny,int index)293 double AdrOp::GetIndexWidth(int ny, int index)
294 {
295 	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
296 	double width=0;
297 	while (ny<0) ny+=uiDimension;
298 	ny=ny%uiDimension;
299 	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
300 	if ((index>0) && (index<(int)uiMax[ny]-1)) width= (dGrid[ny][index+1]-dGrid[ny][index-1])/2*dDeltaUnit;
301 	else if ((index==0) && (index<(int)uiMax[ny]-1)) width=(dGrid[ny][index+1]-dGrid[ny][index])*dDeltaUnit;
302 	else if ((index>0) && (index==(int)uiMax[ny]-1)) width=(dGrid[ny][index]-dGrid[ny][index-1])*dDeltaUnit;
303 	else width= 0;
304 	return width;
305 }
306 
GetIndexCoord(int ny,int index)307 double AdrOp::GetIndexCoord(int ny, int index)
308 {
309 	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
310 	while (ny<0) ny+=uiDimension;
311 	ny=ny%uiDimension;
312 	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
313 	if (index<0) index=0;
314 	if (index>=(int)uiMax[ny]) index=uiMax[ny]-1;
315 	return dGrid[ny][index]*dDeltaUnit;
316 }
317 
GetIndexDelta(int ny,int index)318 double AdrOp::GetIndexDelta(int ny, int index)
319 {
320 	if (index<0) return GetIndexCoord(ny, 0) - GetIndexCoord(ny, 1);
321 	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
322 	if (index>=(int)uiMax[ny]-1) return GetIndexCoord(ny, (int)uiMax[ny]-2) - GetIndexCoord(ny, (int)uiMax[ny]-1);
323 	return GetIndexCoord(ny, index+1) - GetIndexCoord(ny, index);
324 }
325 
326 
Shift(int ny,int step)327 unsigned int AdrOp::Shift(int ny, int step)
328 {
329 	if (bPosSet==false) error->Error(6);
330 	while (ny<0) ny+=uiDimension;
331 	ny=ny%uiDimension;
332 	switch (ny)
333 	{
334 	case 0:
335 		{
336 			iIshift=step;
337 //			if ((int)uiIpos+step<0) iIshift=-2*uiIpos-iIshift;
338 //			else if ((int)uiIpos+step>=(int)uiImax) iIshift=-1*iIshift+2*(uiImax-1-uiIpos);
339 			break;
340 		}
341 	case 1:
342 		{
343 			iJshift=step;
344 //			if ((int)uiJpos+iJshift<0) iJshift=-2*uiJpos-iJshift;
345 //			else if ((int)uiJpos+iJshift>=(int)uiJmax) iJshift=-1*iJshift+2*(uiJmax-1-uiJpos);
346 			break;
347 		}
348 	case 2:
349 		{
350 			iKshift=step;
351 //			if ((int)uiKpos+iKshift<0) iKshift=-2*uiKpos-iKshift;
352 //			else if ((int)uiKpos+iKshift>=(int)uiKmax) iKshift=-1*iKshift+2*(uiKmax-1-uiKpos);
353 			break;
354 		}
355 	}
356 	ADRESSEXPENSE(1,1,0,0,2,3)
357 	return GetPos(iIshift,iJshift,iKshift);
358 }
359 
CheckShift(int ny,int step)360 bool AdrOp::CheckShift(int ny, int step)
361 {
362 	while (ny<0) ny+=uiDimension;
363 	ny=ny%uiDimension;
364 	int shift[3]={0,0,0};
365 	shift[ny]=step;
366 	if (this->CheckPos(uiIpos+shift[0],uiJpos+shift[1],uiKpos+shift[2]))
367 	{
368 		this->Shift(ny,step);
369 		return true;
370 	}
371 	else return false;
372 }
373 
GetShiftedPos(int ny,int step)374 unsigned int AdrOp::GetShiftedPos(int ny, int step)
375 {
376 	if ((ny<0) || (ny>2))
377 		return GetPos(iIshift,iJshift,iKshift);
378 	int pos[3] = {iIshift, iJshift, iKshift};
379 	pos[ny]+=step;
380 	return GetPos(pos[0],pos[1],pos[2]);
381 }
382 
ResetShift()383 void AdrOp::ResetShift()
384 {
385 	iIshift=iJshift=iKshift=0;
386 }
387 
Iterate(int jump)388 unsigned int AdrOp::Iterate(int jump)
389 {
390 	if (abs(jump)>=(int)uiImax) error->Error(4);
391 	i=uiIpos+jump;
392 	if (i>=uiImax)
393 	{
394 		i=i-uiImax;
395 		j=uiJpos+1;
396 		if (j>=uiJmax)
397 		{
398 			j=0;
399 			if (uiDimension==3)
400 			{
401 				k=uiKpos+1;
402 				if (k>=uiKmax) k=0;
403 				uiKpos=k;
404 			}
405 		}
406 		uiIpos=i;
407 		uiJpos=j;
408 		return GetPos();
409 	}
410 	else
411 	{
412 		uiIpos=i;
413 		return GetPos();
414 	}
415 }
416 
GetSize()417 unsigned int AdrOp::GetSize()
418 {
419 	return uiSize;
420 }
421 
422 
SetReflection2Node()423 void AdrOp::SetReflection2Node()
424 {
425 	reflect=true;
426 	uiTypeOffset=0;
427 }
428 
SetReflection2Cell()429 void AdrOp::SetReflection2Cell()
430 {
431 	reflect=true;
432 	uiTypeOffset=1;
433 }
434 
SetReflectionOff()435 void AdrOp::SetReflectionOff()
436 {
437 	reflect=false;
438 }
439 
AddCellAdrOp()440 AdrOp* AdrOp::AddCellAdrOp()
441 {
442 	if (clCellAdr!=NULL) return clCellAdr;
443 	if (uiDimension==3) clCellAdr = new AdrOp(uiImax-1,uiJmax-1,uiKmax-1);
444 	else if (uiDimension==2) clCellAdr = new AdrOp(uiImax-1,uiJmax-1);
445 	else clCellAdr=NULL;
446 	if (clCellAdr!=NULL)
447 	{
448 		clCellAdr->SetPos(0,0,0);
449 		clCellAdr->SetReflection2Cell();
450 	}
451 	iCellShift[0]=iCellShift[1]=iCellShift[2]=0;
452 	return clCellAdr;
453 }
454 
DeleteCellAdrOp()455 AdrOp* AdrOp::DeleteCellAdrOp()
456 {
457 	delete clCellAdr;
458 	clCellAdr=NULL;
459 	return NULL;
460 }
461 
ShiftCell(int ny,int step)462 unsigned int AdrOp::ShiftCell(int ny, int step)
463 {
464 	if (clCellAdr==NULL) error->Error(7);
465 	while (ny<0) ny+=uiDimension;
466 	ny=ny%uiDimension;
467 	iCellShift[ny]=step;
468 	ADRESSEXPENSE(3,1,0,0,1,2)
469 	return clCellAdr->GetPos(uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]);
470 }
471 
ShiftCellCheck(int ny,int step)472 bool AdrOp::ShiftCellCheck(int ny, int step)
473 {
474 	return clCellAdr->CheckShift(ny,step);
475 }
476 
ResetCellShift()477 void AdrOp::ResetCellShift()
478 {
479 	if (clCellAdr==NULL) error->Error(7);
480 	iCellShift[0]=iCellShift[1]=iCellShift[2]=0;
481 }
482 
GetCellPos(bool incShift)483 unsigned int AdrOp::GetCellPos(bool incShift)
484 {
485 	if (bPosSet==false) error->Error(6);
486 	if (clCellAdr==NULL) error->Error(7);
487 #if EXPENSE_LOG==1
488 	if (incShift) ADRESSEXPENSE(3,0,0,0,0,2)
489 #endif
490 		if (incShift) return clCellAdr->GetPos(uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]);
491 		else return clCellAdr->GetPos(uiIpos,uiJpos,uiKpos);
492 }
493 
GetCellPos(int i,int j,int k)494 unsigned int AdrOp::GetCellPos(int i, int j, int k)
495 {
496 	if (bPosSet==false) error->Error(6);
497 	return clCellAdr->GetPos(uiIpos+i,uiJpos+j,uiKpos+k);
498 }
499 
500 
GetShiftCellVolume(int ny,int step)501 double AdrOp::GetShiftCellVolume(int ny, int step)
502 {
503 	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
504 	int uiMax[4]={(int)uiImax-1,(int)uiJmax-1,(int)uiKmax-1,(int)uiLmax-1};
505 	while (ny<0) ny+=uiDimension;
506 	ny=ny%uiDimension;
507 	iCellShift[ny]=step;
508 	int uiPos[4]={(int)uiIpos+iCellShift[0],(int)uiJpos+iCellShift[1],(int)uiKpos+iCellShift[2]};
509 	double dVol=1;
510 	for (unsigned int n=0; n<uiDimension; ++n)
511 	{
512 		if (uiMax[n]>0)
513 		{
514 			while ((uiPos[n]<0) || (uiPos[n]>=uiMax[n]))
515 			{
516 				if (uiPos[n]<0) uiPos[n]=-1*uiPos[n]-1;
517 				else if (uiPos[n]>=uiMax[n]) uiPos[n]=uiMax[n]-(uiPos[n]-uiMax[n]+1);
518 			}
519 			dVol*=(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]])*dDeltaUnit;
520 		}
521 	}
522 	return dVol;
523 }
524 
525 
deltaAdrOp(unsigned int max)526 deltaAdrOp::deltaAdrOp(unsigned int max)
527 {
528 	uiMax=max;
529 }
530 
~deltaAdrOp()531 deltaAdrOp::~deltaAdrOp()
532 {
533 }
534 
SetMax(unsigned int max)535 void deltaAdrOp::SetMax(unsigned int max)
536 {
537 	uiMax=max;
538 }
539 
GetAdr(int pos)540 unsigned int deltaAdrOp::GetAdr(int pos)
541 {
542 	if (uiMax==1) return 0;
543 	if (pos<0)	pos=pos*-1;
544 	else if (pos>(int)uiMax-1) pos=2*(uiMax-1)-pos+1;
545 	if ((pos<0) || (pos>(int)uiMax-1))
546 	{
547 		fprintf(stderr," Error exiting... ");
548 		getchar();
549 		exit(-1);
550 	}
551 	return pos;
552 }
553 
554 
555