#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import cxxfem as fem
import os

fem.debug_level(0)

if __name__ == "__main__":

    # fine mesh
    # pars = {
    #     "Lx": 5.0, # width
    #     "Ly": 1.0, # depth
    #     "Lz": 1.0, # height/thickness
    #     "nx": 60,         
    #     "ny": 30,
    #     "nz": 30
    # }

    # same as python
    pars = {
        "Lx": 5.0, # width
        "Ly": 1.0, # depth
        "Lz": 1.0, # height/thickness
        "nx": 15,         
        "ny": 2,
        "nz": 10
    }

    # mono-element for debug
    # pars = {
    #     "Lx": 5.0, # width
    #     "Ly": 1.0, # depth
    #     "Lz": 1.0, # height/thickness
    #     "nx": 20,         
    #     "ny": 1,
    #     "nz": 1
    # }

    pbl = fem.Problem()
    f = os.path.join(os.path.dirname(__file__), 'parallelepiped.geo')
    pbl.load(f, fem.StringDoubleMap(pars))

    m = fem.Medium(pbl, "Volume", E=100.0, nu=0.3)

    # - mechanical BCs

    # clamping
    bcs = []
    for d in ['X','Y','Z']:
        bcs.append(fem.Dirichlet(pbl, "Left", dof=d, val=0.0))
    
    # Vertical displacement
    # bcs.append(fem.Dirichlet(pbl, "Right", dof='Z', val=-1.0)) # prescribed displ    
    # bcs.append(fem.Neumann(pbl, "Corner", 'Z', -0.1)) # concentrated force (seems OK)
    bcs.append(fem.Neumann(pbl, "Top", 'Z', -1e-2)) # distributed load
    bcs.append(fem.Neumann(pbl, "Top", 'Y', -5e-2)) # distributed load

    solver = fem.Solver(pbl)
    solver.linear_solver = "pardiso"
    solver.solve()

    post = fem.Post(solver)
    # post.build_views(ther=False, smooth=False)
    post.build_views()
    post.show_views(['force_vector','stress_tensor'])
    post.set_camera_3D()
    post.write()    
    post.deform()

    args = fem.parseargs()
    if not args.nogui:    
        post.view()