SCOREC / fep

Finite Element Programming course materials
6 stars 4 forks source link

FEP a4 B matrix #27

Closed startrekman closed 3 years ago

startrekman commented 3 years ago

So I am having some difficulty getting the B matrix formed. So I want to do something like

DenseMatrix B; B.SetSize(dof, idk how to generally get what to put here) For a 4 node quad, I think it would be dofs x 4 using variables from a4_stiffness.cpp CalcDShape(ip, &B)

Where I will use B in B^T D B W_ip to evaluate the integral and evaluate the stiffness matrix. I think I am only getting an error because I do not have the size right.

Thank you for your help

startrekman commented 3 years ago

I tried 4x4 for B for order 1 and that is not working either

startrekman commented 3 years ago

In file included from /gpfs/u/software/erp-spack-install/vDev7a99c49_gcc910_1/linux-centos7-zen/gcc-9. wghx3fztcaz52seh57xzob76tnhqc/include/mfem/mesh/triangle.hpp:16, from /gpfs/u/software/erp-spack-install/vDev7a99c49_gcc910_1/linux-centos7-zen/gcc-9. wghx3fztcaz52seh57xzob76tnhqc/include/mfem/mesh/mesh_headers.hpp:21, from /gpfs/u/software/erp-spack-install/vDev7a99c49_gcc910_1/linux-centos7-zen/gcc-9. wghx3fztcaz52seh57xzob76tnhqc/include/mfem/mfem.hpp:42, from /gpfs/u/software/erp-spack-install/vDev7a99c49_gcc910_1/linux-centos7-zen/gcc-9. wghx3fztcaz52seh57xzob76tnhqc/include/mfem.hpp:2, from a4_element_stiffness.cpp:20: /gpfs/u/software/erp-spack-install/vDev7a99c49_gcc910_1/linux-centos7-zen/gcc-9.1.0/mfem-4.2.0-da2wghx 76tnhqc/include/mfem/fem/fe.hpp:371:41: note: initializing argument 2 of ‘virtual void mfem::FiniteE (const mfem::IntegrationPoint&, mfem::DenseMatrix&) const’ 371 | DenseMatrix &dshape) const = 0; | ~~~^~~~ make: *** [a4_element_stiffness] Error 1

startrekman commented 3 years ago

elmat_diff = 0.; for (int j = 0; j < ir.GetNPoints(); j++) { // loop over quadrature points

const IntegrationPoint& ip = ir.IntPoint(j); et->SetIntPoint(&ip);

// TODO: add the code for computation of element matrix corresponding to the diffusive term // and store it into "elmat_diff" DenseMatrix B; B.SetSize(dofs, 4); FiniteElement::CalcDShape(ip, & B);

}

mortezah commented 3 years ago

@startrekman for a 4 node quad, you will have 4 shape functions. And the derivative will be a 4 by 2 matrix:

  1. 1st column will be the derivatives of the all the shape functions with respect to the 1st parametric confidante, ip.x
  2. 2nd column will be the derivatives of the all the shape functions with respect to the 2nd parametric confidante, ip.y

So in your example above you want to have

DenseMatrix B;
B.SetSize(dofs, 2);
fe->CalcDShape(ip, B);

where fe is the finite element object defined on the element.

startrekman commented 3 years ago

When I try the above code, it still results in an error. The error is as follows:

a4_element_stiffness.cpp: In function ‘int main(int, char**)’: a4_element_stiffness.cpp:109:36: error: cannot call member function ‘virtual void mfem::FiniteElement::CalcDShape(const mfem::IntegrationPoint&, mfem::DenseMatrix&) const’ without object 109 | FiniteElement::CalcDShape(ip, B);

mortezah commented 3 years ago

@startrekman My bad I have updated my comment. You cannot call the class member CalcDShape of the class FiniteElemet using the following line

FiniteElement::CalcDShape(ip, B)

You will need to call it on the object of that class, which is defined a the beginning of the loop in the exercise here https://github.com/SCOREC/fep/blob/d16f0b636671f519bcfb8ea91d07fc04638ad091/a4/a4_element_stiffness.cpp#L77

So the correct way of calling CalcDShape would be

fe->CalcDShape(ip, B);
startrekman commented 3 years ago

Thank you, that seems to be working.