# Gaussian spatial push on 2D surface in 3D geometry Answered

Hi there,

I'm looking for an efficient way to apply mechanical load with gaussian spatial amplitude to avoid edge artifacts for a point source push applied on a 3D anisotropic material (e.g. tendon). I've done this presently by pushing using 5 concentric point sources (cylindrical push) having 5 radii (e.g., 10micron to 50 micron with 10 micron interval). A DO loop was used to apply each push; at each execution the radius of the push cyln enlarges and the amplitude picks a gaussian amplitude reduction factor from a separate file. Since its a surface push applied through DO loop, my method is so slow so that I can't use more resolved amplitude factors to make a smooth push (e.g., by using 50 or 100 radii or circles to mimic a better gaussian scale factor). Anyone has an idea to resolve this?

Thank you

Maju

• Section of code to make the push..

c Import time data of Gaussian
data
hist gauss_data * supergaussian.dat \$push_duration 1.0 \$push_offset
end
c Import spatial data of push
data
hist spatial_data * testGaus.dat \$push_width_scale \$push_amplitude
end

c Push definition
func hist gauss_data 1.0
symb sptsize = 2e-4
symb nsptloop = 5

do push_loop ii 1 \$nsptloop 1 /*\$i\$idx-1 1

symb #get { rxcrd rycrd rzcrd } crdnode \$ii 1
symb #get { push_scale } histval spatial_data \$rxcrd
plod
symb radc = \$ii * \$sptsize / \$nsptloop /* calculate radius of cyln
chek off off
vctr vct\$ii 0. 1. 0.
end
end\$ push_loop

• Hi Maju,

When you say the more resolved simulations are so slow, how slow do you mean? There are a few ways models can be sped up but I would need to take a look at the full code to make suggestions. Are you able to share the files via your preferred cloud sharing service?

Best Regards,

Chloe

• To give little more explanation to the above mentioned problems:
Here below is a 3D simulation with grid size 5 micron (earlier mentioned was 15 micron)

Displacement fields after 333 micro seconds for an isotropic (mu = G = 10k) system is displayed in the below figs those spatially filtered at different bandwidths (previously BW = 500 to 10kHz)

Filtering values are in the title : BW [10 10000] = 10Hz to 10kHz.

As you can see, at lower bandwidth, a square pattern is visible in the diagonal directions. At higher cut off of the low limit, squares were almost invisible but the shape (it should be circle) was lost.

Thank you

Maju

• code used:

`mem mb 4058 50c Model id, title, and basic runtime settingstitl skin 3D modelmp omp 16 *rest noc Variablessymb domain_x = 3e-3symb dermis_thickness = 0.5e-3symb domain_z = 3e-3symb src_x = \$domain_x / 2symb src_z = \$domain_z / 2c Materials - Layersymb density = 1000.0symb dermis_mu = 10000.0symb dermis_G = 1 * \$dermis_muc Pushsymb push_duration = 50E-6 symb push_amplitude = 5000.0symb push_offset = 0E-6 /* \$push_duration / 2.0 * 5.0symb sptsize = 5e-4symb push_width = \$sptsizesymb push_width_scale = \$push_widthc symb cp = 35.0 symb q2 = 0symb dbm = 0 /* delta by musymb q1 = 0.5 * ( \$q2 - \$dbm * \$dermis_mu )symb param_lame_lambda = ( 35.0 * 35.0 - 2.0 ) * \$dermis_mu /* \$density * \$cp * \$cp - 2 * \$shearmod_subq_mu /*( 1500.0 * 1500.0 - 2.0 ) * \$shearmod_subq_musymb c11 = \$param_lame_lambda + 2.0 * \$dermis_mu + \$q2symb c12 = \$param_lame_lambda + \$q1symb c22 = \$param_lame_lambda + 2.0 * \$dermis_musymb c23 = \$param_lame_lambdasymb cp = sqrt ( \$c11 / \$density )symb c44 = \$dermis_musymb c55 = \$dermis_Gsymb c66 = \$dermis_Gc Gridsymb grid_spacing = 5E-6 /* \$est_wavelength / \$elm_per_wavelengthc Outputsymb simtime = 0.5e-3c Meshsymb #keycord x 1 0.0 \$domain_x /* \$water_widthsymb #keycord y 1 0.0 \$dermis_thickness /* \$water_width /*\$skin_thickness \$subq_thickness /*\$halfspace_heightsymb #keycord z 1 0.0 \$domain_zc Get number of indices in x and y for this set of keypointssymb #get { idx } rootmax xsymb #get { jdx } rootmax ysymb #get { kdx } rootmax zc Set grid indices in x and ysymb #keyindx i 1 \$idx 1 \$grid_spacing 1symb #keyindx j 1 \$jdx 1 \$grid_spacing 1symb #keyindx k 1 \$kdx 1 \$grid_spacing 1c Get grid sizesymb nx = \$i\$idxsymb ny = \$j\$jdxsymb nz = \$k\$kdxc Construct grid, enforcing plane strain constraintgrid \$i\$idx \$j\$jdx \$k\$kdx /*pstnc Map grid to geometrygeom keypnt \$idx \$jdx \$kdxaxisform angl defn axismatr cart 0.0 0.0 0.0 0 0 0 endcc Materialcmatrwvsp offtype lean prop dermis \$density * * 0.005stif \$c11 \$c12 \$c12 0 0 0 \$c22\$c23 0 0 0 \$c22 0 00 \$c44 0 0 \$c55 0 \$c66 axis dermis axismatrc hrgl subq_layer_matr BBendc matr prop d0 7000. 160.e9 80.e9siteregn dermis \$i1 \$i2 \$j1 \$j2 \$k1 \$k2endgrphplot matrendc termc Boundary conditionsbounside xmin symmside xmax absrside ymin freeside ymax absrside zmin symmside zmax absr endc Import time data of Gaussiandatac Cloud Scheduler - Replaced this with the line below: hist gauss_data * supergaussian.dat \$push_duration 1.0 \$push_offsethist gauss_data * supergaussian.dat \$push_duration 1.0 \$push_offsetendc Push definitionfunc hist gauss_data 1.0symb #get { ii ij ik rdist } clsnode \$src_x 0 \$src_zsymb iiend = \$ii + 1symb ikend = \$ik + 1plodpdef push_loading func \$push_amplitudevctr vct 0. 1. 0.sdef push_loading vct \$ii \$iiend 1 1 \$ik \$ikendendcalc pres acousticmax aprs none apmxendprcsc proc plot saveexec 1grphplot aprsendend\$ procproc plot 10c termcc Runtime calculationscsymb #get { step } timestep c symb tstep_actual = min ( \$step , \$tstep_max )symb nexec = \$simtime / \$step symb nloops = 100 /* \$simtime / \$restimesymb nexec2 = \$nexec / \$nloopssymb #get { tbeg } wtimedo loopi ii 1 \$nloops 1 /* setting up do loopexec \$nexec2 /* run model for 'nexec2' timestepscc FLXDATO OUTPUTcdatafile out skin_multilayer_results_surface\$iiout1 tmas 1 1 1 1out1 xcrdout1 ycrdout1 zcrdout1 matdout1 matsout1 matpc ---c To export only surface data, append * * 1 1 to these variablesout1 xvel * * 1 1 * * /* * * * * 1 1out1 yvel * * 1 1 * * /* * * * * 1 1 out1 zvel * * 1 1 * * /* * * * * 1 1 out1 xdsp * * 1 1 * * /* * * * * 1 1 out1 ydsp * * 1 1 * * /* * * * * 1 1 out1 zdsp * * 1 1 * * /* * * * * 1 1 endend\$ loopisymb #get { tend } wtimesymb trun = \$tend - \$tbegsymb #get { label } jobnamesymb #save '\$label.symb'stop`
• Hi Maju,

We typically suggest creating point sources by adding a circular primitive with a dummy material and loading the outside surface of that primitive using the material interface. So, for your example you could copy the dermis material e.g. dermis2 which will have the same material properties as dermis. In the site command you would add a point source at the desired location using cyln or sphr (depending if 2D or 3D) and give it material dermis2. You can then load the interface between the materials dermis and dermis2 using the command plod sdf2 which would give you your point source. See the command reference for help on the code syntax.

I hope this helps.

Best Regards,

Chloe

• Thank you for the prompt response. I'll check it out.

Maju

• Hi Chloe,

I've implemented sphere push that resolved the problem of 'slowing down' of my program. Longitudinal wave propagation seemed correct, which shows a nice circular pattern radiating away from the point like source on a 3D geometry with Z axis as depth.

Now, if I wait for shear wave,

then the circular shape of the velocity field (in Z for a push on XY plane in a 3D model) is not preserved. please see the picture taken at 4.24ms. Shear velocity used was 3m/s

• c Push
symb push_duration = 50E-6
symb push_amplitude = 5000.0
symb push_offset = 0E-6 /* \$push_duration / 2.0 * 5.0
symb sptsize = 50e-6
symb push_width = \$sptsize
symb push_width_scale = \$push_width

c Materials - Layer
symb density = 1000.0
symb dermis_mu = 10000.0
symb dermis_G = 1 * \$dermis_mu

symb q2 = 0
symb dbm = 0  /* delta by mu
symb q1 = 0.5 * ( \$q2 - \$dbm * \$dermis_mu )

symb param_lame_lambda = ( 35.0 * 35.0 - 2.0 ) * \$dermis_mu /* \$density * \$cp * \$cp - 2 * \$shearmod_subq_mu  /*( 1500.0 * 1500.0 - 2.0 ) * \$shearmod_subq_mu
symb c11 = \$param_lame_lambda + 2.0 * \$dermis_mu + \$q2
symb c12 = \$param_lame_lambda + \$q1
symb c22 = \$param_lame_lambda + 2.0 * \$dermis_mu
symb c23 = \$param_lame_lambda
symb cp = sqrt ( \$c11 / \$density )

symb c44 = \$dermis_mu
symb c55 = \$dermis_G
symb c66 = \$dermis_G

symb dimxy = 0.015
symb dimz = 0.001
symb #keycord x 1 0 \$dimxy \$dimxy
symb #get { idx } rootmax x

symb #keycord y 1 0 \$dimxy \$dimxy
symb #get { jdx } rootmax y

symb #keycord z 1 0 \$dimz
symb #get { kdx } rootmax z

symb box = 0.00005

symb #keyindx i 1 \$idx 1 \$box 1
symb #keyindx j 1 \$jdx 1 \$box 1
symb #keyindx k 1 \$kdx 1 \$box 1

grid \$i\$idx \$j\$jdx \$k\$kdx

symb domainxy = \$dimxy + \$dimxy
geom
xcrd 0 \$domainxy
ycrd 0 \$domainxy
zcrd 0 \$dimz

matr
wvsp off
type lean
prop dermis \$density * * 0.005
stif \$c11 \$c12 \$c12 0 0 0 \$c22
\$c23 0 0 0 \$c22 0 0
0 \$c44 0 0 \$c55 0 \$c66
c c    hrgl subq_layer_matr BB
end
c matr
c     wvsp on
c     type elas
c     prop steel 1000. 110 3
c     end
matr
wvsp on
type elas
prop steel2 \$density \$cp 3
end

site
regn dermis
sphr steel2 \$dimxy  \$dimxy 0 0.0003
end

c Boundary conditions
boun
side xmin absr
side xmax absr
side ymin absr
side ymax absr
side zmin free
side zmax absr
end
grph
plot matr
end

term
data
hist gauss_data * supergaussian.dat \$push_duration 1.0 \$push_offset
end
func hist gauss_data 1.0

plod
pdef pld1 func \$push_amplitude
spot spot1 \$dimxy \$dimxy 0
sdf2 pld1 spot1 dermis steel2
end

calc
pres acoustic
max aprs none apmx
disp
end

prcs

proc plot save

exec 1

grph
plot aprs
end

end\$ proc

proc plot 5

term

symb simtime = 0.005
c
c Runtime calculations
c
symb #get { step } timestep

symb nexec = \$simtime / \$step
symb nloops = 100
symb nexec2 = \$nexec / \$nloops

symb #get { tbeg } wtime

do loopi ii 1 \$nloops 1 /* setting up do loop
exec \$nexec2 /* run model for 'nexec2' timesteps
c
c FLXDATO OUTPUT
c
data
file out skin_multilayer_results_surface\$ii
out1 tmas 1 1 1 1
out1 xcrd
out1 ycrd
out1 zcrd
out1 matd
out1 mats
out1 matp
c ---
c To export only surface data, append  * * 1 1  to these variables
out1 xvel * * * * 1 1
out1 yvel * * * * 1 1
out1 zvel * * * * 1 1
out1 xdsp * * * * 1 1
out1 ydsp * * * * 1 1
out1 zdsp * * * * 1 1
c    out1 xcrl * * * * * *
c    out1 ycrl * * * * * *
c    out1 zcrl * * * * * *
end
end\$ loopi

symb #get { tend } wtime
symb trun = \$tend - \$tbeg

symb #get { label } jobname
symb #save '\$label.symb'

stop

• Hi Maju,

I think the absorbing boundaries may be introducing artefacts at the corners. To combat this you can bend your model into a circle so the absorbing boundary is a smooth surface. This will also ensure the point source surface is smooth (instead of stepped due to the mesh). Here is an example that shows you the syntax for the bend command.

Best Regards,

Chloe