Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81411466
region_intersect.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Sep 5, 22:21
Size
11 KB
Mime Type
text/x-c
Expires
Sat, Sep 7, 22:21 (2 d)
Engine
blob
Format
Raw Data
Handle
20562668
Attached To
rLAMMPS lammps
region_intersect.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <stdlib.h>
#include <string.h>
#include "region_intersect.h"
#include "domain.h"
#include "error.h"
#include "force.h"
using
namespace
LAMMPS_NS
;
/* ---------------------------------------------------------------------- */
RegIntersect
::
RegIntersect
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Region
(
lmp
,
narg
,
arg
)
{
if
(
narg
<
5
)
error
->
all
(
FLERR
,
"Illegal region command"
);
int
n
=
force
->
inumeric
(
FLERR
,
arg
[
2
]);
if
(
n
<
2
)
error
->
all
(
FLERR
,
"Illegal region command"
);
options
(
narg
-
(
n
+
3
),
&
arg
[
n
+
3
]);
// build list of regions to intersect
// store sub-region IDs in idsub
idsub
=
new
char
*
[
n
];
list
=
new
int
[
n
];
nregion
=
0
;
int
m
,
iregion
;
for
(
int
iarg
=
0
;
iarg
<
n
;
iarg
++
)
{
m
=
strlen
(
arg
[
iarg
+
3
])
+
1
;
idsub
[
nregion
]
=
new
char
[
m
];
strcpy
(
idsub
[
nregion
],
arg
[
iarg
+
3
]);
iregion
=
domain
->
find_region
(
idsub
[
nregion
]);
if
(
iregion
==
-
1
)
error
->
all
(
FLERR
,
"Region intersect region ID does not exist"
);
list
[
nregion
++
]
=
iregion
;
}
// this region is variable shape or dynamic if any of sub-regions are
Region
**
regions
=
domain
->
regions
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
{
if
(
regions
[
list
[
ilist
]]
->
varshape
)
varshape
=
1
;
if
(
regions
[
list
[
ilist
]]
->
dynamic
)
dynamic
=
1
;
}
// extent of intersection of regions
// has bounding box if interior and any sub-region has bounding box
bboxflag
=
0
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
if
(
regions
[
list
[
ilist
]]
->
bboxflag
==
1
)
bboxflag
=
1
;
if
(
!
interior
)
bboxflag
=
0
;
if
(
bboxflag
)
{
int
first
=
1
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
{
if
(
regions
[
list
[
ilist
]]
->
bboxflag
==
0
)
continue
;
if
(
first
)
{
extent_xlo
=
regions
[
list
[
ilist
]]
->
extent_xlo
;
extent_ylo
=
regions
[
list
[
ilist
]]
->
extent_ylo
;
extent_zlo
=
regions
[
list
[
ilist
]]
->
extent_zlo
;
extent_xhi
=
regions
[
list
[
ilist
]]
->
extent_xhi
;
extent_yhi
=
regions
[
list
[
ilist
]]
->
extent_yhi
;
extent_zhi
=
regions
[
list
[
ilist
]]
->
extent_zhi
;
first
=
0
;
}
extent_xlo
=
MAX
(
extent_xlo
,
regions
[
list
[
ilist
]]
->
extent_xlo
);
extent_ylo
=
MAX
(
extent_ylo
,
regions
[
list
[
ilist
]]
->
extent_ylo
);
extent_zlo
=
MAX
(
extent_zlo
,
regions
[
list
[
ilist
]]
->
extent_zlo
);
extent_xhi
=
MIN
(
extent_xhi
,
regions
[
list
[
ilist
]]
->
extent_xhi
);
extent_yhi
=
MIN
(
extent_yhi
,
regions
[
list
[
ilist
]]
->
extent_yhi
);
extent_zhi
=
MIN
(
extent_zhi
,
regions
[
list
[
ilist
]]
->
extent_zhi
);
}
}
// possible contacts = sum of possible contacts in all sub-regions
// for near contacts and touching contacts
cmax
=
0
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
cmax
+=
regions
[
list
[
ilist
]]
->
cmax
;
contact
=
new
Contact
[
cmax
];
tmax
=
0
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
{
if
(
interior
)
tmax
+=
regions
[
list
[
ilist
]]
->
tmax
;
else
tmax
++
;
}
}
/* ---------------------------------------------------------------------- */
RegIntersect
::~
RegIntersect
()
{
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
delete
[]
idsub
[
ilist
];
delete
[]
idsub
;
delete
[]
list
;
delete
[]
contact
;
}
/* ---------------------------------------------------------------------- */
void
RegIntersect
::
init
()
{
Region
::
init
();
// re-build list of sub-regions in case other regions were deleted
// error if a sub-region was deleted
int
iregion
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
{
iregion
=
domain
->
find_region
(
idsub
[
ilist
]);
if
(
iregion
==
-
1
)
error
->
all
(
FLERR
,
"Region union region ID does not exist"
);
list
[
ilist
]
=
iregion
;
}
// init the sub-regions
Region
**
regions
=
domain
->
regions
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
regions
[
list
[
ilist
]]
->
init
();
}
/* ----------------------------------------------------------------------
inside = 1 if x,y,z is match() with all sub-regions
else inside = 0
------------------------------------------------------------------------- */
int
RegIntersect
::
inside
(
double
x
,
double
y
,
double
z
)
{
int
ilist
;
Region
**
regions
=
domain
->
regions
;
for
(
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
if
(
!
regions
[
list
[
ilist
]]
->
match
(
x
,
y
,
z
))
break
;
if
(
ilist
==
nregion
)
return
1
;
return
0
;
}
/* ----------------------------------------------------------------------
compute contacts with interior of intersection of sub-regions
(1) compute contacts in each sub-region
(2) only keep a contact if surface point is match() to all other regions
------------------------------------------------------------------------- */
int
RegIntersect
::
surface_interior
(
double
*
x
,
double
cutoff
)
{
int
m
,
ilist
,
jlist
,
iregion
,
jregion
,
ncontacts
;
double
xs
,
ys
,
zs
;
Region
**
regions
=
domain
->
regions
;
int
n
=
0
;
int
walloffset
=
0
;
for
(
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
{
iregion
=
list
[
ilist
];
ncontacts
=
regions
[
iregion
]
->
surface
(
x
[
0
],
x
[
1
],
x
[
2
],
cutoff
);
for
(
m
=
0
;
m
<
ncontacts
;
m
++
)
{
xs
=
x
[
0
]
-
regions
[
iregion
]
->
contact
[
m
].
delx
;
ys
=
x
[
1
]
-
regions
[
iregion
]
->
contact
[
m
].
dely
;
zs
=
x
[
2
]
-
regions
[
iregion
]
->
contact
[
m
].
delz
;
for
(
jlist
=
0
;
jlist
<
nregion
;
jlist
++
)
{
if
(
jlist
==
ilist
)
continue
;
jregion
=
list
[
jlist
];
if
(
!
regions
[
jregion
]
->
match
(
xs
,
ys
,
zs
))
break
;
}
if
(
jlist
==
nregion
)
{
contact
[
n
].
r
=
regions
[
iregion
]
->
contact
[
m
].
r
;
contact
[
n
].
radius
=
regions
[
iregion
]
->
contact
[
m
].
radius
;
contact
[
n
].
delx
=
regions
[
iregion
]
->
contact
[
m
].
delx
;
contact
[
n
].
dely
=
regions
[
iregion
]
->
contact
[
m
].
dely
;
contact
[
n
].
delz
=
regions
[
iregion
]
->
contact
[
m
].
delz
;
contact
[
n
].
iwall
=
regions
[
iregion
]
->
contact
[
m
].
iwall
+
walloffset
;
contact
[
n
].
varflag
=
regions
[
iregion
]
->
contact
[
m
].
varflag
;
n
++
;
}
}
// increment by cmax instead of tmax to insure
// possible wall IDs for sub-regions are non overlapping
walloffset
+=
regions
[
iregion
]
->
cmax
;
}
return
n
;
}
/* ----------------------------------------------------------------------
compute contacts with interior of intersection of sub-regions
(1) flip interior/exterior flag of each sub-region
(2) compute contacts in each sub-region
(3) only keep a contact if surface point is not match() to all other regions
(4) flip interior/exterior flags back to original settings
this is effectively same algorithm as surface_interior() for RegUnion
------------------------------------------------------------------------- */
int
RegIntersect
::
surface_exterior
(
double
*
x
,
double
cutoff
)
{
int
m
,
ilist
,
jlist
,
iregion
,
jregion
,
ncontacts
;
double
xs
,
ys
,
zs
;
Region
**
regions
=
domain
->
regions
;
int
n
=
0
;
for
(
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
regions
[
list
[
ilist
]]
->
interior
^=
1
;
for
(
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
{
iregion
=
list
[
ilist
];
ncontacts
=
regions
[
iregion
]
->
surface
(
x
[
0
],
x
[
1
],
x
[
2
],
cutoff
);
for
(
m
=
0
;
m
<
ncontacts
;
m
++
)
{
xs
=
x
[
0
]
-
regions
[
iregion
]
->
contact
[
m
].
delx
;
ys
=
x
[
1
]
-
regions
[
iregion
]
->
contact
[
m
].
dely
;
zs
=
x
[
2
]
-
regions
[
iregion
]
->
contact
[
m
].
delz
;
for
(
jlist
=
0
;
jlist
<
nregion
;
jlist
++
)
{
if
(
jlist
==
ilist
)
continue
;
jregion
=
list
[
jlist
];
if
(
regions
[
jregion
]
->
match
(
xs
,
ys
,
zs
))
break
;
}
if
(
jlist
==
nregion
)
{
contact
[
n
].
r
=
regions
[
iregion
]
->
contact
[
m
].
r
;
contact
[
n
].
radius
=
regions
[
iregion
]
->
contact
[
m
].
radius
;
contact
[
n
].
delx
=
regions
[
iregion
]
->
contact
[
m
].
delx
;
contact
[
n
].
dely
=
regions
[
iregion
]
->
contact
[
m
].
dely
;
contact
[
n
].
delz
=
regions
[
iregion
]
->
contact
[
m
].
delz
;
contact
[
n
].
iwall
=
ilist
;
contact
[
n
].
varflag
=
regions
[
iregion
]
->
contact
[
m
].
varflag
;
n
++
;
}
}
}
for
(
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
regions
[
list
[
ilist
]]
->
interior
^=
1
;
return
n
;
}
/* ----------------------------------------------------------------------
change region shape of all sub-regions
------------------------------------------------------------------------- */
void
RegIntersect
::
shape_update
()
{
Region
**
regions
=
domain
->
regions
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
regions
[
list
[
ilist
]]
->
shape_update
();
}
/* ----------------------------------------------------------------------
move/rotate all sub-regions
------------------------------------------------------------------------- */
void
RegIntersect
::
pretransform
()
{
Region
**
regions
=
domain
->
regions
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
regions
[
list
[
ilist
]]
->
pretransform
();
}
/* ----------------------------------------------------------------------
get translational/angular velocities of all subregions
------------------------------------------------------------------------- */
void
RegIntersect
::
set_velocity
()
{
Region
**
regions
=
domain
->
regions
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
regions
[
list
[
ilist
]]
->
set_velocity
();
}
/* ----------------------------------------------------------------------
increment length of restart buffer based on region info
used by restart of fix/wall/gran/region
------------------------------------------------------------------------- */
void
RegIntersect
::
length_restart_string
(
int
&
n
)
{
n
+=
sizeof
(
int
)
+
strlen
(
id
)
+
1
+
sizeof
(
int
)
+
strlen
(
style
)
+
1
+
sizeof
(
int
);
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
domain
->
regions
[
list
[
ilist
]]
->
length_restart_string
(
n
);
}
/* ----------------------------------------------------------------------
region writes its current position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
void
RegIntersect
::
write_restart
(
FILE
*
fp
)
{
int
sizeid
=
(
strlen
(
id
)
+
1
);
int
sizestyle
=
(
strlen
(
style
)
+
1
);
fwrite
(
&
sizeid
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
id
,
1
,
sizeid
,
fp
);
fwrite
(
&
sizestyle
,
sizeof
(
int
),
1
,
fp
);
fwrite
(
style
,
1
,
sizestyle
,
fp
);
fwrite
(
&
nregion
,
sizeof
(
int
),
1
,
fp
);
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
){
domain
->
regions
[
list
[
ilist
]]
->
write_restart
(
fp
);
}
}
/* ----------------------------------------------------------------------
region reads its previous position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
int
RegIntersect
::
restart
(
char
*
buf
,
int
&
n
)
{
int
sizeid
=
buf
[
n
];
n
+=
sizeof
(
int
);
char
*
restart_id
=
new
char
[
sizeid
];
for
(
int
i
=
0
;
i
<
sizeid
;
i
++
)
restart_id
[
i
]
=
buf
[
n
++
];
if
(
strcmp
(
restart_id
,
id
)
!=
0
)
return
0
;
int
sizestyle
=
buf
[
n
];
n
+=
sizeof
(
int
);
char
*
restart_style
=
new
char
[
sizestyle
];
for
(
int
i
=
0
;
i
<
sizestyle
;
i
++
)
restart_style
[
i
]
=
buf
[
n
++
];
if
(
strcmp
(
restart_style
,
style
)
!=
0
)
return
0
;
int
restart_nreg
=
buf
[
n
];
n
+=
sizeof
(
int
);
if
(
restart_nreg
!=
nregion
)
return
0
;
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
){
if
(
!
domain
->
regions
[
list
[
ilist
]]
->
restart
(
buf
,
n
)){
return
0
;
}
}
return
1
;
}
/* ----------------------------------------------------------------------
set prev vector to zero
------------------------------------------------------------------------- */
void
RegIntersect
::
reset_vel
()
{
for
(
int
ilist
=
0
;
ilist
<
nregion
;
ilist
++
)
domain
->
regions
[
list
[
ilist
]]
->
reset_vel
();
}
Event Timeline
Log In to Comment