code: def bernstein_space(order, nsd):
      if nsd > 3:
          raise RuntimeError(
              "Bernstein only implemented in 1D, 2D, and 3D"
          )
      sum = 0
      basis = []
      coeff = []
  
      if nsd == 1:
          b1, b2 = x, 1 - x
          for o1 in range(0, order + 1):
              for o2 in range(0, order + 1):
                  if o1 + o2 == order:
                      aij = Symbol("a_%d_%d" % (o1, o2))
                      sum += (
                          aij
                          * binomial(order, o1)
                          * pow(b1, o1)
                          * pow(b2, o2)
                      )
                      basis.append(
                          binomial(order, o1)
                          * pow(b1, o1)
                          * pow(b2, o2)
                      )
                      coeff.append(aij)
  
      if nsd == 2:
          b1, b2, b3 = x, y, 1 - x - y
          for o1 in range(0, order + 1):
              for o2 in range(0, order + 1):
                  for o3 in range(0, order + 1):
                      if o1 + o2 + o3 == order:
                          aij = Symbol("a_%d_%d_%d" % (o1, o2, o3))
                          fac = factorial(order) / (
                              factorial(o1)
                              * factorial(o2)
                              * factorial(o3)
                          )
                          sum += (
                              aij
                              * fac
                              * pow(b1, o1)
                              * pow(b2, o2)
                              * pow(b3, o3)
                          )
                          basis.append(
                              fac
                              * pow(b1, o1)
                              * pow(b2, o2)
                              * pow(b3, o3)
                          )
                          coeff.append(aij)
  
      if nsd == 3:
          b1, b2, b3, b4 = x, y, z, 1 - x - y - z
          for o1 in range(0, order + 1):
              for o2 in range(0, order + 1):
                  for o3 in range(0, order + 1):
                      for o4 in range(0, order + 1):
                          if o1 + o2 + o3 + o4 == order:
                              aij = Symbol(
                                  "a_%d_%d_%d_%d" % (o1, o2, o3, o4)
                              )
                              fac = factorial(order) / (
                                  factorial(o1)
                                  * factorial(o2)
                                  * factorial(o3)
                                  * factorial(o4)
                              )
                              sum += (
                                  aij
                                  * fac
                                  * pow(b1, o1)
                                  * pow(b2, o2)
                                  * pow(b3, o3)
                                  * pow(b4, o4)
                              )
                              basis.append(
                                  fac
                                  * pow(b1, o1)
                                  * pow(b2, o2)
                                  * pow(b3, o3)
                                  * pow(b4, o4)
                              )
                              coeff.append(aij)
  
      return sum, coeff, basis
  
file_path: sympy/examples/advanced/fem.py
imports: from sympy import symbols, Symbol, factorial, Rational, zeros, eye, integrate, diff, pprint, reduced
label_desc: Use image node_id 1 for calling a global function with example usage: bernstein_space(order, nsd) and returns: sum,  coeff,  basis