How to build stiffness matrix in two dimensions?

[thehalfspace]

I am trying to assemble stiffness matrix as a sparse matrix in two dimensions, similar to what is shown in one dimension.

I know the current problem calculates K.u at a local level for efficiency as shown here, but in quasi-static problem the matrix solver is a very slow step.

I am trying to do something like this:

# lagrange derivatives and weights as obtained from GetGLL
H = H[:]
Ht = Ht[:]
wgll2 = wgll*wgll' 
jac = dx_dxi*dy_deta   # jacobian in 2D

kglob = zeros(nglob,nglob)

# Ke is elemental stiffness
Ke = zeros(NGLL*NGLL,NGLL*NGLL)

for i in each element in x and y:
    n = 0  # iterator
    for xe in 1:NGLL:         # loop for each spectral node within one element
        for ye in 1:NGLL:
             n = n+1
              for ng in 1:NGLL*NGLL:    # loop for shape functions
                      Ke[n,ng] += mu * wgll2[xe,ye] * jac * H[ng] * Ht[ng]

       end for ng, ye, xe

      ig = iglob[:,:,i][:]
      kglob[ig,ig] += Ke
end for i

This doesn't give me correct answer. I am confused due to the symmetry of shape functions in x and y.

[jpampuero]

You may find a better explanation of the operations here and (for irregular grids) in section 3.6 of my lecture notes.

[thehalfspace]

Thanks. The lecture notes were super helpful. I am able to assemble the stiffness as sparse matrix now.

[jpampuero]

I am glad it helped. Please consider contributing your code to semlab, if you think others can benefit from it. If you make it available somewhere else, let me know.

[thehalfspace]

I'll definitely contribute to SEMLAB once I am done testing my stuff. The multiplication (Ksparse*d) is much faster than assembling the local contributions using loop.

I have done the sparse matrix assembly here. I'll try to rewrite it in matlab and submit a pr for SEMLAB.

I am working on the earthquake cycle code which I will make available on Github once it is fully functional. Here is the prototype for now.

[jpampuero]

Great work in Julia. The most efficient versions we have achieved in SEMLAB so far are: sem2d_NewmarkA1_scec2_vectorized.m sem2d_NewmarkA1_scec2_vectorized_parfor.m We have used those optimizations to model cycles too. It will be interesting to see if a sparse matrix implementation performs better in matlab. Performance comparisons may depend on the language.