made adjustments to support potential new boundary handling scheme; not sure if I'll go through with it so I may roll these changes back

This commit is contained in:
Paul M. Sutter 2025-05-22 18:36:01 -04:00
parent 23d665f7bd
commit 05030b24e8
2 changed files with 122 additions and 35 deletions

View file

@ -133,7 +133,7 @@ int main(int argc, char *argv[]) {
DL c[d] = ((float)b[d]+0.5)*width;
printf("Counting particles in subdomain...\n");
/* Assign temporary array*/
nvpbuf = 0; /* Number of particles to tesselate, includng ibuffer */
nvpbuf = 0; /* Number of particles to tesselate, includng buffer and guards */
nvp = 0; /* Without the buffer */
for (i=0; i<np; i++) {
isitinbuf = 1;
@ -237,49 +237,121 @@ int main(int argc, char *argv[]) {
/* Add guard points */
printf("Adding guard points to subdomain faces...\n");
// first the internal boundaries between subdomains
printf("Adding guard points to internal subdomain faces...\n");
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
for (j=0; j<NGUARD+1; j++) {
/* Bottom */
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = -width2 - g;
nvpall++;
if (b[2] != 0) { // we will do the outer edges next
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = -width2 - g;
nvpall++;
}
/* Top */
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = width2 + g;
nvpall++;
}
if (b[2] != numdiv-1) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = width2 + g;
nvpall++;
}
}
}
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 - g;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
for (j=0; j<NGUARD+1; j++) {
if (b[1] != 0) { // we will do the outer edges next
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 - g;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = width2 + g;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
if (b[1] != numdiv-1) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = width2 + g;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
}
}
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 - g;
parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
for (j=0; j<NGUARD+1; j++) {
if (b[0] != 0) { // we will do the outer edges next
parts[3*nvpall] = -width2 - g;
parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
parts[3*nvpall] = width2 + g;
parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
if (b[0] != numdiv-1) {
parts[3*nvpall] = width2 + g;
parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
}
}
/*
// now the external faces of the domain
// (i'm tired and can't think of a better way to do this)
isObs = 0;
if (isObs != 1) {
printf("Adding guard points to external domain faces...\n");
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
/* Bottom */
if (b[2] == 0) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = -width2 - g;
nvpall++;
}
/* Top */
if (b[2] == numdiv-1) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = width2 + g;
nvpall++;
}
}
}
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
for (j=0; j<NGUARD+1; j++) {
if (b[1] == 0) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 - g;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
if (b[1] == numdiv-1) {
parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = width2 + g;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
}
}
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
if (b[0] == 0) {
parts[3*nvpall] = -width2 - g;
parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
if (b[0] == numdiv-1) {
parts[3*nvpall] = width2 + g;
parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++;
}
}
}
} // end if isObs != -1
// TEST TEMP
testFile = fopen("test.txt", "a");
for (i=0; i<nvpall; i++){
@ -293,7 +365,6 @@ int main(int argc, char *argv[]) {
parts[3*i+2]);
}
fclose(testFile);
*/
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=nvpbuf;i<nvpall;i++) {

View file

@ -130,6 +130,21 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
log.write("Done!\n")
log.close()
# add tracers along the bounding box to contain the tessellation
log.write("\nAdding guards along box faces...\n")
nGuard = 42
# load in file
# store number of particles
# store number of guards
# place guards
# re-save file
if jobSuccessful(logFile, "Done!\n"):
endTime = time.time()
@ -497,6 +512,7 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
nadjPerPart = np.fromfile(File, dtype=np.int32,count=numPart)
for p in range(numPart):
nin = np.fromfile(File, dtype=np.int32, count=1)[0]
part[p].nadjs = nin
for n in range(nin):
pAdj = np.fromfile(File, dtype=np.int32, count=1)[0]
if (p < pAdj):