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  

9 comments

  • 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
    pdef push_loading$ii func $push_scale
    vctr vct$ii 0. 1. 0.
    cyln push_loading$ii vct$ii stnd y 0 $grid_spacing 0.0 0.0 $radc
    end
    end$ push_loop

    Comment actions Permalink
  • 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

    Comment actions Permalink
  • 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

    Comment actions Permalink
  • code used:



    mem mb 4058 50

    c Model id, title, and basic runtime settings

    titl skin 3D model
    mp omp 16 *
    rest no

    c Variables

    symb domain_x = 3e-3
    symb dermis_thickness = 0.5e-3
    symb domain_z = 3e-3
    symb src_x = $domain_x / 2
    symb src_z = $domain_z / 2
    c Materials - Layer
    symb density = 1000.0
    symb dermis_mu = 10000.0
    symb dermis_G = 1 * $dermis_mu
    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 = 5e-4
    symb push_width = $sptsize
    symb push_width_scale = $push_width

    c symb cp = 35.0
    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


    c Grid


    symb grid_spacing = 5E-6 /* $est_wavelength / $elm_per_wavelength


    c Output

    symb simtime = 0.5e-3




    c Mesh

    symb #keycord x 1 0.0 $domain_x /* $water_width
    symb #keycord y 1 0.0 $dermis_thickness /* $water_width /*$skin_thickness $subq_thickness /*$halfspace_height
    symb #keycord z 1 0.0 $domain_z

    c Get number of indices in x and y for this set of keypoints
    symb #get { idx } rootmax x
    symb #get { jdx } rootmax y
    symb #get { kdx } rootmax z

    c Set grid indices in x and y
    symb #keyindx i 1 $idx 1 $grid_spacing 1
    symb #keyindx j 1 $jdx 1 $grid_spacing 1
    symb #keyindx k 1 $kdx 1 $grid_spacing 1

    c Get grid size
    symb nx = $i$idx
    symb ny = $j$jdx
    symb nz = $k$kdx

    c Construct grid, enforcing plane strain constraint
    grid $i$idx $j$jdx $k$kdx /*pstn

    c Map grid to geometry
    geom keypnt $idx $jdx $kdx

    axis
    form angl
    defn axismatr cart 0.0 0.0 0.0 0 0 0
    end


    c
    c Material
    c

    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
    axis dermis axismatr
    c hrgl subq_layer_matr BB
    end
    c matr prop d0 7000. 160.e9 80.e9

    site
    regn dermis $i1 $i2 $j1 $j2 $k1 $k2
    end

    grph
    plot matr
    end
    c term

    c Boundary conditions
    boun
    side xmin symm
    side xmax absr
    side ymin free
    side ymax absr
    side zmin symm
    side zmax absr
    end


    c Import time data of Gaussian
    data

    c Cloud Scheduler - Replaced this with the line below: hist gauss_data * supergaussian.dat $push_duration 1.0 $push_offset
    hist gauss_data * supergaussian.dat $push_duration 1.0 $push_offset

    end


    c Push definition
    func hist gauss_data 1.0
    symb #get { ii ij ik rdist } clsnode $src_x 0 $src_z
    symb iiend = $ii + 1
    symb ikend = $ik + 1
    plod
    pdef push_loading func $push_amplitude
    vctr vct 0. 1. 0.
    sdef push_loading vct $ii $iiend 1 1 $ik $ikend
    end







    calc
    pres acoustic
    max aprs none apmx
    end
    prcs
    c
    proc plot save

    exec 1

    grph

    plot aprs
    end
    end$ proc

    proc plot 10
    c
    term





    c
    c Runtime calculations
    c
    symb #get { step } timestep
    c symb tstep_actual = min ( $step , $tstep_max )
    symb nexec = $simtime / $step
    symb nloops = 100 /* $simtime / $restime
    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 * * /* * * * * 1 1
    out1 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

    end








    end$ loopi

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


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

    stop
    Comment actions Permalink
  • 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

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

    Maju

    Comment actions Permalink
  • 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

    Comment actions Permalink
  • 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

    Comment actions Permalink
  • 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

    Comment actions Permalink

Please sign in to leave a comment.

Didn't find what you were looking for?

New post