reset('')
forget()
load("my_classes.sage")
load("functions.sage")

# choose basis functions
var('xhat yhat zhat');
assume(xhat>-1,xhat<1,yhat>-1,yhat<1)
phi = get_triangular_basis(xhat,yhat)
#phi = get_quadrilateral_basis(xhat,yhat)

# define integration limits
ndof = 4;
limY1 = -1; limY2 = 1
limX1 = -1; limX2 = 1
if (len(phi) == 3):
    ndof = 3;
    limY1 = 0; limY2 = 1-xhat;
    limX1 = 0; limX2 = 1;


# build nodal arrays
var('x1 x2 x3 x4')
var('y1 y2 y3 y4')
XX1 = SVECT(x1, y1, 0);XX2 = SVECT(x2, y2, 0);XX3 = SVECT(x3, y3, 0);XX4 = SVECT(x4, y4, 0);
p  = [XX1, XX2, XX3, XX4];

# calculate 2d element jacobian and cartesian mapped shape function gradients
Adjac2d,dphi = getAnalytic2dDjacsDphi(p,phi);
Adjac2d = Adjac2d.maxima_methods().rootscontract().full_simplify().canonicalize_radical().collect_common_factors();

# arbitrary functions
dof_map = [0,1,2,3];
var('f_avg fx_avg fy_avg g_avg h_avg eu ev code_djac')
var('f1 f2 f3 f4')
fi = [f1, f2, f3, f4];
f = expand_function(fi,phi,dof_map);

# area using Jacobian
Area = (integral( integral( Adjac2d, yhat,limY1,limY2), xhat,limX1,limX2)).maxima_methods().rootscontract().full_simplify().canonicalize_radical().collect_common_factors()

print "\nAdjac2d: ",Adjac2d
print "Area: ",Area

###################################################################
# integrations
###################################################################
result = [0]*len(phi);

###################################################################
print "\n**integral(constant * phi_i)**"
con = 3/Area
print "double con = one_3 * area"
for j in range(len(phi)):
    result[j]  = integral( integral( phi[j] * Adjac2d * con, yhat,limY1,limY2), xhat,limX1,limX2)
    print "integral[",j,"] = con * (",result[j].maxima_methods().rootscontract().full_simplify().canonicalize_radical().collect_common_factors(),")"

###################################################################
print "\n**lumped mass matrices - integral(phi_i * (f_i * phi_i))**"
con = 6/Area
print "double con = one_6 * area"
for j in range(len(phi)):
    result[j]  = integral( integral( phi[j] * (fi[j] * phi[j]) * Adjac2d * con, yhat,limY1,limY2), xhat,limX1,limX2)
    result[j] = result[j].maxima_methods().rootscontract().full_simplify().canonicalize_radical().collect_common_factors()
    print "integral[",j,"] = con * (",result[j],")"

###################################################################
print "\n** what I think you guys are integrating **"
con = 3/Area
print "double con = one_3 * area"
for j in range(len(phi)):
    result[j]  = integral( integral( phi[j] * (fi[j]) * Adjac2d * con, yhat,limY1,limY2), xhat,limX1,limX2)
    result[j] = result[j].maxima_methods().rootscontract().full_simplify().canonicalize_radical().collect_common_factors()
    print "integral[",j,"] = con * (",result[j],")"